diff --git a/.gitmodules b/.gitmodules index d098a48b7f..e4e3f9c097 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,7 +1,7 @@ [submodule "physics/Radiation/RRTMGP/rte-rrtmgp"] path = physics/Radiation/RRTMGP/rte-rrtmgp url = https://github.com/NCAR/rte-rrtmgp - branch = main + branch = production/GFS.v17 [submodule "physics/MP/TEMPO/TEMPO"] path = physics/MP/TEMPO/TEMPO url = https://github.com/NCAR/TEMPO diff --git a/physics/MP/GFDL/v1_2019/gfdl_cloud_microphys_mod.F90 b/physics/MP/GFDL/v1_2019/gfdl_cloud_microphys_mod.F90 index c29a116fbc..ad0304074e 100644 --- a/physics/MP/GFDL/v1_2019/gfdl_cloud_microphys_mod.F90 +++ b/physics/MP/GFDL/v1_2019/gfdl_cloud_microphys_mod.F90 @@ -2837,7 +2837,7 @@ subroutine lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, q, precip, m1, mono) k0 = ktop do k = ktop, kbot - do n = k0, kbot + n_loop: do n = k0, kbot if (ze (k) <= zt (n) .and. ze (k) >= zt (n + 1)) then pl = (zt (n) - ze (k)) / dz (n) if (zt (n + 1) <= ze (k + 1)) then @@ -2847,7 +2847,7 @@ subroutine lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, q, precip, m1, mono) a4 (4, n) * r3 * (pr * (pr + pl) + pl ** 2) qm (k) = qm (k) * (ze (k) - ze (k + 1)) k0 = n - goto 555 + exit n_loop else qm (k) = (ze (k) - zt (n + 1)) * (a4 (2, n) + 0.5 * (a4 (4, n) + & a4 (3, n) - a4 (2, n)) * (1. + pl) - a4 (4, n) * (r3 * (1. + pl * (1. + pl)))) @@ -2862,15 +2862,14 @@ subroutine lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, q, precip, m1, mono) qm (k) = qm (k) + delz * (a4 (2, m) + 0.5 * esl * & (a4 (3, m) - a4 (2, m) + a4 (4, m) * (1. - r23 * esl))) k0 = m - goto 555 + exit n_loop endif enddo endif - goto 555 + exit n_loop endif endif - enddo - 555 continue + enddo n_loop enddo m1 (ktop) = q (ktop) - qm (ktop) diff --git a/physics/MP/GFDL/v3_2022/gfdl_cloud_microphys_v3_mod.F90 b/physics/MP/GFDL/v3_2022/gfdl_cloud_microphys_v3_mod.F90 index b15f2efd98..55aaa9e0a2 100644 --- a/physics/MP/GFDL/v3_2022/gfdl_cloud_microphys_v3_mod.F90 +++ b/physics/MP/GFDL/v3_2022/gfdl_cloud_microphys_v3_mod.F90 @@ -4542,7 +4542,7 @@ subroutine lagrangian_fall (ks, ke, zs, ze, zt, dp, q, precip, m1) k0 = ks do k = ks, ke - do n = k0, ke + n_loop: do n = k0, ke if (ze (k) .le. zt (n) .and. ze (k) .ge. zt (n + 1)) then pl = (zt (n) - ze (k)) / dz (n) if (zt (n + 1) .le. ze (k + 1)) then @@ -4552,7 +4552,7 @@ subroutine lagrangian_fall (ks, ke, zs, ze, zt, dp, q, precip, m1) a4 (4, n) * r3 * (pr * (pr + pl) + pl ** 2) qm (k) = qm (k) * (ze (k) - ze (k + 1)) k0 = n - goto 555 + exit n_loop else qm (k) = (ze (k) - zt (n + 1)) * (a4 (2, n) + 0.5 * (a4 (4, n) + & a4 (3, n) - a4 (2, n)) * (1. + pl) - a4 (4, n) * (r3 * (1. + pl * (1. + pl)))) @@ -4567,15 +4567,14 @@ subroutine lagrangian_fall (ks, ke, zs, ze, zt, dp, q, precip, m1) qm (k) = qm (k) + delz * (a4 (2, m) + 0.5 * esl * & (a4 (3, m) - a4 (2, m) + a4 (4, m) * (1. - r23 * esl))) k0 = m - goto 555 + exit n_loop endif enddo endif - goto 555 + exit n_loop endif endif - enddo - 555 continue + enddo n_loop enddo m1 (ks) = q (ks) - qm (ks) diff --git a/physics/MP/Morrison_Gettelman/aer_cloud.F b/physics/MP/Morrison_Gettelman/aer_cloud.F index 36bdf47acd..f38a2e7566 100644 --- a/physics/MP/Morrison_Gettelman/aer_cloud.F +++ b/physics/MP/Morrison_Gettelman/aer_cloud.F @@ -1543,7 +1543,7 @@ subroutine activate (wparc,ndroplet,smax,nmodes, ! ! *** perform bisection ************************************************* ! -20 do 30 i=1,maxit_par + do i=1,maxit_par x3 = 0.5*(x1+x2) ! if (ntot .gt. zero_par) then @@ -1562,15 +1562,15 @@ subroutine activate (wparc,ndroplet,smax,nmodes, x1 = x3 endif ! - if (abs(x2-x1) .le. eps_par*x1) goto 40 + if (abs(x2-x1) .le. eps_par*x1) exit niter = i -30 continue + end do ! ! *** converged ; return ************************************************ ! -40 x3 = 0.5*(x1+x2) + x3 = 0.5*(x1+x2) ! if (ntot .gt. zero_par) then call sintegral (x2,ndrpl,sinteg1,sinteg2,wparcel,nmodes, @@ -1826,25 +1826,29 @@ subroutine gauleg (x,w,n) m=(n+1)/2d0 xm=0.5d0*(x2+x1) xl=0.5d0*(x2-x1) - do 12 i=1,m + + do i=1,m z=cos(pi_par*(i-.25d0)/(n+.5d0)) - 1 continue + do p1=1.d0 p2=0.d0 - do 11 j=1,n - p3=p2 - p2=p1 - p1=((2.d0*j-1.)*z*p2-(j-1.d0)*p3)/j - 11 continue + do j=1,n + p3=p2 + p2=p1 + p1=((2.d0*j-1.)*z*p2-(j-1.d0)*p3)/j + end do + pp=n*(z*p1-p2)/(z*z-1.d0) z1=z z=z1-p1/pp - if(abs(z-z1).gt.eps_par)go to 1 + if(abs(z-z1).gt.eps_par) exit + end do + x(i)=xm-xl*z x(n+1-i)=xm+xl*z w(i)=2.d0*xl/((1.d0-z*z)*pp*pp) w(n+1-i)=w(i) - 12 continue + end do end subroutine gauleg !C======================================================================= @@ -2922,28 +2926,22 @@ real*8 function cubicint_ice(y, y1, y2, a, b) real*8 :: A_, B_, a0, a1, a2, a3, d, AUX if (y .le. y1) then - d=a - goto 5065 - end if - - if (y .ge. y2) then - d=b - goto 5065 - end if - - - AUX=y2-y1 - A_=6d0*(a-b)/(AUX*AUX*AUX) - B_=a+(A_*(y1*y1*y1)/6d0)-(A_*(y1*y1)*y2*0.5d0) - - a0=B_ - a1=A_*y1*y2 - a2=-A_*(y1+y2)*0.5d0 - a3=A_/3d0 - d=a0+(a1*y)+(a2*y*y)+(a3*y*y*y) - + d=a + elseif (y .ge. y2) then + d=b + else + AUX=y2-y1 + A_=6d0*(a-b)/(AUX*AUX*AUX) + B_=a+(A_*(y1*y1*y1)/6d0)-(A_*(y1*y1)*y2*0.5d0) + + a0=B_ + a1=A_*y1*y2 + a2=-A_*(y1+y2)*0.5d0 + a3=A_/3d0 + d=a0+(a1*y)+(a2*y*y)+(a3*y*y*y) + endif - 5065 cubicint_ice=d + cubicint_ice=d end function cubicint_ice @@ -2958,26 +2956,20 @@ real*8 function dcubicint_ice(y, y1, y2, a, b) real*8 :: A_, a0, a1, a2, a3, d, AUX if (y .le. y1) then - d=0 - goto 5065 - end if - - if (y .ge. y2) then - d=0 - goto 5065 - end if - - - AUX=y2-y1 - A_=6d0*(a-b)/(AUX*AUX*AUX) - - a1=A_*y1*y2 - a2=-A_*(y1+y2)*0.5d0 - a3=A_/3d0 - d=(a1)+(2d0*a2*y)+(3d0*a3*y*y) + d=0 + elseif (y .ge. y2) then + d=0 + else + AUX=y2-y1 + A_=6d0*(a-b)/(AUX*AUX*AUX) + a1=A_*y1*y2 + a2=-A_*(y1+y2)*0.5d0 + a3=A_/3d0 + d=(a1)+(2d0*a2*y)+(3d0*a3*y*y) + endif - 5065 dcubicint_ice=d + dcubicint_ice=d end function dcubicint_ice diff --git a/physics/MP/Morrison_Gettelman/cldwat2m_micro.F b/physics/MP/Morrison_Gettelman/cldwat2m_micro.F index a80790eb6e..d0e59bedec 100644 --- a/physics/MP/Morrison_Gettelman/cldwat2m_micro.F +++ b/physics/MP/Morrison_Gettelman/cldwat2m_micro.F @@ -4745,11 +4745,10 @@ subroutine findsp1 (lchnk, ncol, q, t, p, tsp, qsp) end do if (dtm < dttol .and. dqm < dqtol) then - go to 10 + exit endif end do - 10 continue error_found = .false. if (dtm > dttol .or. dqm > dqtol) then @@ -5023,11 +5022,10 @@ subroutine findsp1_water (lchnk, ncol, q, t, p, tsp, qsp) end do if (dtm < dttol .and. dqm < dqtol) then - go to 10 + exit endif end do - 10 continue error_found = .false. if (dtm > dttol .or. dqm > dqtol) then @@ -5439,7 +5437,8 @@ FUNCTION GAMMA(X) Y = Y + ONE ELSE RES=XINF - GOTO 900 + GAMMA = RES + RETURN ENDIF ENDIF !---------------------------------------------------------------------- @@ -5453,7 +5452,8 @@ FUNCTION GAMMA(X) RES = ONE/Y ELSE RES = XINF - GOTO 900 + GAMMA = RES + RETURN ENDIF ELSEIF(Y.LT.TWELVE)THEN Y1 = Y @@ -5510,7 +5510,8 @@ FUNCTION GAMMA(X) RES = EXP(SUM) ELSE RES = XINF - GOTO 900 + GAMMA = RES + RETURN ENDIF ENDIF !---------------------------------------------------------------------- @@ -5518,7 +5519,7 @@ FUNCTION GAMMA(X) !---------------------------------------------------------------------- IF(PARITY) RES = -RES IF(FACT.NE.ONE) RES = FACT/RES -900 GAMMA = RES + GAMMA = RES !D900 DGAMMA = RES RETURN ! ---------- LAST LINE OF GAMMA ---------- diff --git a/physics/MP/Morrison_Gettelman/wv_saturation.F b/physics/MP/Morrison_Gettelman/wv_saturation.F index b12b76a917..6ec4919a66 100644 --- a/physics/MP/Morrison_Gettelman/wv_saturation.F +++ b/physics/MP/Morrison_Gettelman/wv_saturation.F @@ -211,7 +211,7 @@ subroutine gestbl(tmn ,tmx ,trice ,ip ,epsil , latvap ,latice , & return ! -9000 format('GESTBL: FATAL ERROR ********************************* & +9000 format('GESTBL: FATAL ERROR ********************************* & &',/, ' TMAX AND TMIN REQUIRE A LARGER DIMENSION ON THE LENGTH', & & ' OF THE SATURATION VAPOR PRESSURE TABLE ESTBL(PLENEST)',/, & & ' TMAX, TMIN, AND PLENEST => ', 2f7.2, i3) @@ -407,11 +407,38 @@ subroutine aqsatd(t, p, es, qs, gam, ii, ilen, kk, kstart, kend) ! accurate to within 1 percent for 173.16 < t < 373.16 ! trinv = 0.0_kp - if ((.not. icephs) .or. (ttrice == 0.0_kp)) go to 10 - trinv = 1.0_kp/ttrice + if ((.not. icephs) .or. (ttrice == 0.0_kp)) then ! - do k=kstart,kend - do i=1,ilen +! No icephs or water to ice transition +! + do k=kstart,kend + do i=1,ilen +! +! Account for change of hlatv with t above freezing where +! constant slope is given by -2369 j/(kg c) = cpv - cw +! + hlatvp = hlatv - 2369.0_kp*(t(i,k)-tmelt) + if (icephs) then + hlatsb = hlatv + hlatf + else + hlatsb = hlatv + end if + if (t(i,k) < tmelt) then + hltalt = hlatsb + else + hltalt = hlatvp + end if + desdt = hltalt*es(i,k)/(rgasv*t(i,k)*t(i,k)) + gam(i,k) = hltalt*qs(i,k)*p(i,k)*desdt/(cp*es(i,k)*(p(i,k) & + & - omeps*es(i,k))) + if (qs(i,k) == 1.0_kp) gam(i,k) = 0.0_kp + end do + end do + else + trinv = 1.0_kp/ttrice +! + do k=kstart,kend + do i=1,ilen ! ! Weighting of hlat accounts for transition from water to ice ! polynomial expression approximates difference between es over @@ -420,58 +447,32 @@ subroutine aqsatd(t, p, es, qs, gam, ii, ilen, kk, kstart, kend) ! range from ice to water also accounting for change of hlatv with t ! above freezing where constant slope is given by -2369 j/(kg c) =cpv - cw ! - tc = t(i,k) - tmelt - lflg = (tc >= -ttrice .and. tc < 0.0_kp) - weight = min(-tc*trinv,1.0_kp) - hlatsb = hlatv + weight*hlatf - hlatvp = hlatv - 2369.0_kp*tc - if (t(i,k) < tmelt) then - hltalt = hlatsb - else - hltalt = hlatvp - end if - if (lflg) then - tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3) + tc*(pcf(4) & - & + tc*pcf(5)))) - else - tterm = 0.0_kp - end if - desdt = hltalt*es(i,k)/(rgasv*t(i,k)*t(i,k)) + tterm*trinv - gam(i,k) = hltalt*qs(i,k)*p(i,k)*desdt/(cp*es(i,k)*(p(i,k) & - & - omeps*es(i,k))) - if (qs(i,k) == 1.0_kp) gam(i,k) = 0.0_kp - end do - end do -! - go to 20 -! -! No icephs or water to ice transition -! -10 do k=kstart,kend - do i=1,ilen -! -! Account for change of hlatv with t above freezing where -! constant slope is given by -2369 j/(kg c) = cpv - cw -! - hlatvp = hlatv - 2369.0_kp*(t(i,k)-tmelt) - if (icephs) then - hlatsb = hlatv + hlatf - else - hlatsb = hlatv - end if - if (t(i,k) < tmelt) then - hltalt = hlatsb - else - hltalt = hlatvp - end if - desdt = hltalt*es(i,k)/(rgasv*t(i,k)*t(i,k)) - gam(i,k) = hltalt*qs(i,k)*p(i,k)*desdt/(cp*es(i,k)*(p(i,k) & - & - omeps*es(i,k))) - if (qs(i,k) == 1.0_kp) gam(i,k) = 0.0_kp + tc = t(i,k) - tmelt + lflg = (tc >= -ttrice .and. tc < 0.0_kp) + weight = min(-tc*trinv,1.0_kp) + hlatsb = hlatv + weight*hlatf + hlatvp = hlatv - 2369.0_kp*tc + if (t(i,k) < tmelt) then + hltalt = hlatsb + else + hltalt = hlatvp + end if + if (lflg) then + tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3) + tc*(pcf(4) + & + tc*pcf(5)))) + else + tterm = 0.0_kp + end if + desdt = hltalt*es(i,k)/(rgasv*t(i,k)*t(i,k)) + & + tterm*trinv + gam(i,k) = hltalt*qs(i,k)*p(i,k)*desdt/ + & (cp*es(i,k)*(p(i,k) - omeps*es(i,k))) + if (qs(i,k) == 1.0_kp) gam(i,k) = 0.0_kp + end do end do - end do + end if ! -20 return + return end subroutine aqsatd !>\ingroup wv_saturation_mod @@ -539,9 +540,34 @@ subroutine vqsatd(t ,p ,es ,qs ,gam , len ) ! accurate to within 1 percent for 173.16 < t < 373.16 ! trinv = 0.0_kp - if ((.not. icephs) .or. (ttrice.eq.0.0_kp)) go to 10 - trinv = 1.0_kp/ttrice - do i=1,len + if ((.not. icephs) .or. (ttrice.eq.0.0_kp)) then +! +! No icephs or water to ice transition +! + do i=1,len +! +! Account for change of hlatv with t above freezing where +! constant slope is given by -2369 j/(kg c) = cpv - cw +! + hlatvp = hlatv - 2369.0_kp*(t(i)-tmelt) + if (icephs) then + hlatsb = hlatv + hlatf + else + hlatsb = hlatv + end if + if (t(i) < tmelt) then + hltalt = hlatsb + else + hltalt = hlatvp + end if + desdt = hltalt*es(i)/(rgasv*t(i)*t(i)) + gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) + & - omeps*es(i))) + if (qs(i) == 1.0_kp) gam(i) = 0.0_kp + end do + else + trinv = 1.0_kp/ttrice + do i=1,len ! ! Weighting of hlat accounts for transition from water to ice ! polynomial expression approximates difference between es over @@ -550,50 +576,28 @@ subroutine vqsatd(t ,p ,es ,qs ,gam , len ) ! range from ice to water also accounting for change of hlatv with t ! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw ! - tc = t(i) - tmelt - lflg = (tc >= -ttrice .and. tc < 0.0_kp) - weight = min(-tc*trinv,1.0_kp) - hlatsb = hlatv + weight*hlatf - hlatvp = hlatv - 2369.0_kp*tc - if (t(i) < tmelt) then - hltalt = hlatsb - else - hltalt = hlatvp - end if - if (lflg) then - tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3) + tc*(pcf(4) & + tc = t(i) - tmelt + lflg = (tc >= -ttrice .and. tc < 0.0_kp) + weight = min(-tc*trinv,1.0_kp) + hlatsb = hlatv + weight*hlatf + hlatvp = hlatv - 2369.0_kp*tc + if (t(i) < tmelt) then + hltalt = hlatsb + else + hltalt = hlatvp + end if + if (lflg) then + tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3) + tc*(pcf(4) & + tc*pcf(5)))) - else - tterm = 0.0_kp - end if - desdt = hltalt*es(i)/(rgasv*t(i)*t(i)) + tterm*trinv - gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) - omeps*es(i))) - if (qs(i) == 1.0_kp) gam(i) = 0.0_kp - end do - return -! -! No icephs or water to ice transition -! -10 do i=1,len -! -! Account for change of hlatv with t above freezing where -! constant slope is given by -2369 j/(kg c) = cpv - cw -! - hlatvp = hlatv - 2369.0_kp*(t(i)-tmelt) - if (icephs) then - hlatsb = hlatv + hlatf - else - hlatsb = hlatv - end if - if (t(i) < tmelt) then - hltalt = hlatsb - else - hltalt = hlatvp + else + tterm = 0.0_kp + end if + desdt = hltalt*es(i)/(rgasv*t(i)*t(i)) + tterm*trinv + gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) + & - omeps*es(i))) + if (qs(i) == 1.0_kp) gam(i) = 0.0_kp + end do end if - desdt = hltalt*es(i)/(rgasv*t(i)*t(i)) - gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) - omeps*es(i))) - if (qs(i) == 1.0_kp) gam(i) = 0.0_kp - end do ! return ! @@ -1098,9 +1102,37 @@ subroutine vqsatd2(t ,p ,es ,qs ,dqsdt , len ) ! accurate to within 1 percent for 173.16 < t < 373.16 ! trinv = 0.0_kp - if ((.not. icephs) .or. (ttrice == 0.0_kp)) go to 10 - trinv = 1.0_kp/ttrice - do i=1,len + if ((.not. icephs) .or. (ttrice == 0.0_kp)) then +! +! No icephs or water to ice transition +! + do i=1,len +! +! Account for change of hlatv with t above freezing where +! constant slope is given by -2369 j/(kg c) = cpv - cw +! + hlatvp = hlatv - 2369.0_kp*(t(i)-tmelt) + if (icephs) then + hlatsb = hlatv + hlatf + else + hlatsb = hlatv + end if + if (t(i) < tmelt) then + hltalt = hlatsb + else + hltalt = hlatvp + end if + desdt = hltalt*es(i)/(rgasv*t(i)*t(i)) + gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) & + & - omeps*es(i))) + if (qs(i) == 1.0_kp) gam(i) = 0.0_kp + + dqsdt(i) = (cp/hltalt)*gam(i) + + end do + else + trinv = 1.0_kp/ttrice + do i=1,len ! ! Weighting of hlat accounts for transition from water to ice ! polynomial expression approximates difference between es over @@ -1109,57 +1141,31 @@ subroutine vqsatd2(t ,p ,es ,qs ,dqsdt , len ) ! range from ice to water also accounting for change of hlatv with t ! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw ! - tc = t(i) - tmelt - lflg = (tc >= -ttrice .and. tc < 0.0_kp) - weight = min(-tc*trinv,1.0_kp) - hlatsb = hlatv + weight*hlatf - hlatvp = hlatv - 2369.0_kp*tc - if (t(i) < tmelt) then - hltalt = hlatsb - else - hltalt = hlatvp - end if - if (lflg) then - tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3) + tc*(pcf(4) & + tc = t(i) - tmelt + lflg = (tc >= -ttrice .and. tc < 0.0_kp) + weight = min(-tc*trinv,1.0_kp) + hlatsb = hlatv + weight*hlatf + hlatvp = hlatv - 2369.0_kp*tc + if (t(i) < tmelt) then + hltalt = hlatsb + else + hltalt = hlatvp + end if + if (lflg) then + tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3) + tc*(pcf(4) & & + tc*pcf(5)))) - else - tterm = 0.0_kp - end if - desdt = hltalt*es(i)/(rgasv*t(i)*t(i)) + tterm*trinv - gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) - omeps*es(i))) - if (qs(i) == 1.0_kp) gam(i) = 0.0_kp + else + tterm = 0.0_kp + end if + desdt = hltalt*es(i)/(rgasv*t(i)*t(i)) + tterm*trinv + gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) & + & - omeps*es(i))) + if (qs(i) == 1.0_kp) gam(i) = 0.0_kp - dqsdt(i) = (cp/hltalt)*gam(i) + dqsdt(i) = (cp/hltalt)*gam(i) - end do - return -! -! No icephs or water to ice transition -! -10 do i=1,len -! -! Account for change of hlatv with t above freezing where -! constant slope is given by -2369 j/(kg c) = cpv - cw -! - hlatvp = hlatv - 2369.0_kp*(t(i)-tmelt) - if (icephs) then - hlatsb = hlatv + hlatf - else - hlatsb = hlatv - end if - if (t(i) < tmelt) then - hltalt = hlatsb - else - hltalt = hlatvp + end do end if - desdt = hltalt*es(i)/(rgasv*t(i)*t(i)) - gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) - omeps*es(i))) - if (qs(i) == 1.0_kp) gam(i) = 0.0_kp - - dqsdt(i) = (cp/hltalt)*gam(i) - - end do -! return ! end subroutine vqsatd2 @@ -1247,8 +1253,35 @@ subroutine vqsatd2_single(t ,p ,es ,qs ,dqsdt) ! accurate to within 1 percent for 173.16 < t < 373.16 ! trinv = 0.0_kp - if ((.not. icephs) .or. (ttrice == 0.0_kp)) go to 10 - trinv = 1.0_kp/ttrice + if ((.not. icephs) .or. (ttrice == 0.0_kp)) then +! +! No icephs or water to ice transition +! +!10 do i=1,len +! +! Account for change of hlatv with t above freezing where +! constant slope is given by -2369 j/(kg c) = cpv - cw +! + hlatvp = hlatv - 2369.0_kp*(t-tmelt) + if (icephs) then + hlatsb = hlatv + hlatf + else + hlatsb = hlatv + end if + if (t < tmelt) then + hltalt = hlatsb + else + hltalt = hlatvp + end if + desdt = hltalt*es/(rgasv*t*t) + gam = hltalt*qs*p*desdt/(cp*es*(p - omeps*es)) + if (qs == 1.0_kp) gam = 0.0_kp + + dqsdt = (cp/hltalt)*gam + +! end do + else + trinv = 1.0_kp/ttrice ! do i=1,len ! @@ -1259,61 +1292,31 @@ subroutine vqsatd2_single(t ,p ,es ,qs ,dqsdt) ! range from ice to water also accounting for change of hlatv with t ! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw ! - tc = t - tmelt - lflg = (tc >= -ttrice .and. tc < 0.0_kp) - weight = min(-tc*trinv,1.0_kp) - hlatsb = hlatv + weight*hlatf - hlatvp = hlatv - 2369.0_kp*tc - if (t < tmelt) then - hltalt = hlatsb - else - hltalt = hlatvp - end if - if (lflg) then - tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3) + tc*(pcf(4) & + tc = t - tmelt + lflg = (tc >= -ttrice .and. tc < 0.0_kp) + weight = min(-tc*trinv,1.0_kp) + hlatsb = hlatv + weight*hlatf + hlatvp = hlatv - 2369.0_kp*tc + if (t < tmelt) then + hltalt = hlatsb + else + hltalt = hlatvp + end if + if (lflg) then + tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3) + tc*(pcf(4) & & + tc*pcf(5)))) - else - tterm = 0.0_kp - end if - desdt = hltalt*es/(rgasv*t*t) + tterm*trinv - gam = hltalt*qs*p*desdt/(cp*es*(p - omeps*es)) - if (qs == 1.0_kp) gam = 0.0_kp + else + tterm = 0.0_kp + end if + desdt = hltalt*es/(rgasv*t*t) + tterm*trinv + gam = hltalt*qs*p*desdt/(cp*es*(p - omeps*es)) + if (qs == 1.0_kp) gam = 0.0_kp - dqsdt = (cp/hltalt)*gam + dqsdt = (cp/hltalt)*gam ! end do - return -! -! No icephs or water to ice transition -! - -10 continue - -!10 do i=1,len -! -! Account for change of hlatv with t above freezing where -! constant slope is given by -2369 j/(kg c) = cpv - cw -! - hlatvp = hlatv - 2369.0_kp*(t-tmelt) - if (icephs) then - hlatsb = hlatv + hlatf - else - hlatsb = hlatv end if - if (t < tmelt) then - hltalt = hlatsb - else - hltalt = hlatvp - end if - desdt = hltalt*es/(rgasv*t*t) - gam = hltalt*qs*p*desdt/(cp*es*(p - omeps*es)) - if (qs == 1.0_kp) gam = 0.0_kp - - dqsdt = (cp/hltalt)*gam - - -! end do -! + return ! end subroutine vqsatd2_single @@ -1405,42 +1408,43 @@ subroutine gffgch(t ,es ,tmelt ,itype ) end if ! - if(t < (tmelt - tr) .and. itype == 1) go to 10 + if (.not. (t < (tmelt - tr) .and. itype == 1)) then ! ! Water ! - ps = 1013.246_kp - ts = 373.16_kp - e1 = 11.344_kp*(1.0_kp - t/ts) - e2 = -3.49149_kp*(ts/t - 1.0_kp) - f1 = -7.90298_kp*(ts/t - 1.0_kp) - f2 = 5.02808_kp*log10(ts/t) - f3 = -1.3816_kp*(10.0_kp**e1 - 1.0_kp)/10000000.0_kp - f4 = 8.1328_kp*(10.0_kp**e2 - 1.0_kp)/1000.0_kp - f5 = log10(ps) - f = f1 + f2 + f3 + f4 + f5 - es = (10.0_kp**f)*100.0_kp - eswtr = es + ps = 1013.246_kp + ts = 373.16_kp + e1 = 11.344_kp*(1.0_kp - t/ts) + e2 = -3.49149_kp*(ts/t - 1.0_kp) + f1 = -7.90298_kp*(ts/t - 1.0_kp) + f2 = 5.02808_kp*log10(ts/t) + f3 = -1.3816_kp*(10.0_kp**e1 - 1.0_kp)/10000000.0_kp + f4 = 8.1328_kp*(10.0_kp**e2 - 1.0_kp)/1000.0_kp + f5 = log10(ps) + f = f1 + f2 + f3 + f4 + f5 + es = (10.0_kp**f)*100.0_kp + eswtr = es + end if ! - if(t >= tmelt .or. itype == 0) go to 20 + if (.not. (t >= tmelt .or. itype == 0)) then ! ! Ice ! -10 continue - t0 = tmelt - term1 = 2.01889049_kp/(t0/t) - term2 = 3.56654_kp*log(t0/t) - term3 = 20.947031_kp*(t0/t) - es = 575.185606e10_kp*exp(-(term1 + term2 + term3)) + t0 = tmelt + term1 = 2.01889049_kp/(t0/t) + term2 = 3.56654_kp*log(t0/t) + term3 = 20.947031_kp*(t0/t) + es = 575.185606e10_kp*exp(-(term1 + term2 + term3)) ! - if (t < (tmelt - tr)) go to 20 + if (.not. (t < (tmelt - tr))) then ! ! Weighted transition between water and ice ! - weight = min((tmelt - t)/tr,1.0_kp) - es = weight*es + (1.0_kp - weight)*eswtr + weight = min((tmelt - t)/tr,1.0_kp) + es = weight*es + (1.0_kp - weight)*eswtr + end if + end if ! -20 continue itype = itypo return ! @@ -1549,7 +1553,7 @@ subroutine vqsatd2_ice_single(t ,p ,es ,qs ,dqsdt) ! end subroutine vqsatd2_ice_single -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! end module wv_saturation -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! diff --git a/physics/Radiation/RRTMGP/rte-rrtmgp b/physics/Radiation/RRTMGP/rte-rrtmgp index 41c5fcd950..5564070c1e 160000 --- a/physics/Radiation/RRTMGP/rte-rrtmgp +++ b/physics/Radiation/RRTMGP/rte-rrtmgp @@ -1 +1 @@ -Subproject commit 41c5fcd950fed09b8afe186dede266824eca7fd3 +Subproject commit 5564070c1eb1b9e4102e74d7e255b80f8324a3b8 diff --git a/physics/SFC_Models/Land/RUC/module_sf_ruclsm.F90 b/physics/SFC_Models/Land/RUC/module_sf_ruclsm.F90 index aa4f34d1fa..9b6a05717e 100644 --- a/physics/SFC_Models/Land/RUC/module_sf_ruclsm.F90 +++ b/physics/SFC_Models/Land/RUC/module_sf_ruclsm.F90 @@ -795,7 +795,7 @@ SUBROUTINE LSMRUC(xlat,xlon, & do k=2,nzs if(zsmain(k).ge.0.4_kind_phys) then NROOT=K - goto 111 + exit endif enddo ELSE @@ -810,11 +810,10 @@ SUBROUTINE LSMRUC(xlat,xlon, & do k=2,nzs if(zsmain(k).ge.1.1_kind_phys) then NROOT=K - goto 111 + exit endif enddo ENDIF - 111 continue !----- IF (debug_print ) THEN @@ -1523,10 +1522,10 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j, & !--- input varia if(snhei.gt.0.0081_kind_phys*rhowater/rhosn) then !*** Update snow density for current temperature (Koren et al 1999,doi:10.1029/1999JD900232.) BSN=delt/3600._kind_phys*c1sn*exp(0.08_kind_phys*min(zero,tsnav)-c2sn*rhosn*1.e-3_kind_phys) - if(bsn*snwe*100._kind_phys.lt.1.e-4_kind_phys) goto 777 - XSN=rhosn*(exp(bsn*snwe*100._kind_phys)-one)/(bsn*snwe*100._kind_phys) - rhosn=MIN(MAX(58.8_kind_phys,XSN),500._kind_phys) - 777 continue + if(bsn*snwe*100._kind_phys.ge.1.e-4_kind_phys) then + XSN=rhosn*(exp(bsn*snwe*100._kind_phys)-one)/(bsn*snwe*100._kind_phys) + rhosn=MIN(MAX(58.8_kind_phys,XSN),500._kind_phys) + endif endif !-- snow_mosaic from the previous time step @@ -2322,13 +2321,15 @@ FUNCTION QSN(TN,T) R=(TN-173.15_kind_dbl_prec)/.05_kind_dbl_prec+one I=INT(R) - IF(I.GE.1) goto 10 - I=1 - R=1. - 10 IF(I.LE.5000) GOTO 20 - I=5000 - R=5001._kind_dbl_prec - 20 R1=T(I) + + if (I .LT. 1) then + I = 1 + R = 1._kind_dbl_prec + else if (I .GT. 5000) then + I = 5000 + R = 5001._kind_dbl_prec + end if + R1=T(I) R2=R-I QSN=(T(I+1)-R1)*R2 + R1 !----------------------------------------------------------------------- @@ -4857,7 +4858,7 @@ SUBROUTINE SOILTEMP( debug_print,xlat,xlon,testptlat,testptlon,& !****************************************************************************** cotso(1)=zero rhtso(1)=TSO(NZS) - DO 33 K=1,NZS2 + DO K=1,NZS2 KN=NZS-K K1=2*KN-3 X1=DTDZS(K1)*THDIF(KN-1) @@ -4867,7 +4868,7 @@ SUBROUTINE SOILTEMP( debug_print,xlat,xlon,testptlat,testptlon,& DENOM=1.+X1+X2-X2*cotso(K) cotso(K+1)=X1/DENOM rhtso(K+1)=(FT+X2*rhtso(K))/DENOM - 33 CONTINUE + END DO !************************************************************************ !--- THE HEAT BALANCE EQUATION (Smirnova et al., 1996, EQ. 21,26)