diff --git a/physics/GFS_surface_composites.F90 b/physics/GFS_surface_composites.F90 index cd5f3db11..9636eb384 100644 --- a/physics/GFS_surface_composites.F90 +++ b/physics/GFS_surface_composites.F90 @@ -123,7 +123,7 @@ subroutine GFS_surface_composites_pre_run (im, frac_grid, flag_cice, cplflx, lan if (cice(i) < one) then wet(i) = .true. ! tsfco(i) = tgice - tsfco(i) = max(tisfc(i), tgice) + if (.not. cplflx) tsfco(i) = max(tisfc(i), tgice) ! tsfco(i) = max((tsfc(i) - cice(i)*tisfc(i)) & ! / (one - cice(i)), tgice) endif diff --git a/physics/drag_suite.F90 b/physics/drag_suite.F90 index 269bf0b3a..080bee156 100644 --- a/physics/drag_suite.F90 +++ b/physics/drag_suite.F90 @@ -332,11 +332,11 @@ subroutine drag_suite_run( & & hpbl(im), & & slmsk(im) real(kind=kind_phys), dimension(im) :: govrth,xland - real(kind=kind_phys), dimension(im,km) :: dz2 + !real(kind=kind_phys), dimension(im,km) :: dz2 real(kind=kind_phys) :: tauwavex0,tauwavey0, & & XNBV,density,tvcon,hpbl2 integer :: kpbl2,kvar - real(kind=kind_phys), dimension(im,km+1) :: zq ! = PHII/g + !real(kind=kind_phys), dimension(im,km+1) :: zq ! = PHII/g real(kind=kind_phys), dimension(im,km) :: zl ! = PHIL/g !SPP @@ -413,10 +413,10 @@ subroutine drag_suite_run( & ! local variables ! integer :: i,j,k,lcap,lcapp1,nwd,idir, & - klcap,kp1,ikount,kk + klcap,kp1 ! - real(kind=kind_phys) :: rcs,rclcs,csg,fdir,cleff,cleff_ss,cs, & - rcsks,wdir,ti,rdz,temp,tem2,dw2,shr2, & + real(kind=kind_phys) :: rcs,csg,fdir,cleff,cleff_ss,cs, & + rcsks,wdir,ti,rdz,tem2,dw2,shr2, & bvf2,rdelks,wtkbj,tem,gfobnv,hd,fro, & rim,temc,tem1,efact,temv,dtaux,dtauy, & dtauxb,dtauyb,eng0,eng1 @@ -442,7 +442,6 @@ subroutine drag_suite_run( & coefm(im),coefm_ss(im) ! integer :: kbl(im),klowtop(im) - logical :: iope integer,parameter :: mdir=8 !integer :: nwdir(mdir) !data nwdir/6,7,5,8,2,3,1,4/ @@ -661,6 +660,17 @@ subroutine drag_suite_run( & enddo enddo ! +! calculate mid-layer height (zl), interface height (zq), and layer depth (dz2). +! + !zq=0. + do k = kts,km + do i = its,im + !zq(i,k+1) = PHII(i,k+1)*g_inv + !dz2(i,k) = (PHII(i,k+1)-PHII(i,k))*g_inv + zl(i,k) = PHIL(i,k)*g_inv + enddo + enddo +! ! determine reference level: maximum of 2*var and pbl heights ! do i = its,im @@ -895,7 +905,6 @@ subroutine drag_suite_run( & density=1.2 utendwave=0. vtendwave=0. - zq=0. ! IF ( (gwd_opt_ss .EQ. 1).and.(ss_taper.GT.1.E-02) ) THEN if (me==master) print *,"in Drag Suite: Running small-scale gravity wave drag" @@ -914,14 +923,6 @@ subroutine drag_suite_run( & thvx(i,k) = thx(i,k)*tvcon enddo enddo - ! Calculate mid-layer height (zl), interface height (zq), and layer depth (dz2). - do k = kts,km - do i = its,im - zq(i,k+1) = PHII(i,k+1)*g_inv - dz2(i,k) = (PHII(i,k+1)-PHII(i,k))*g_inv - zl(i,k) = PHIL(i,k)*g_inv - enddo - enddo do i=its,im hpbl2 = hpbl(i)+10. @@ -1027,19 +1028,6 @@ subroutine drag_suite_run( & utendform=0. vtendform=0. - zq=0. - - IF ( (gwd_opt_ss .NE. 1).and.(ss_taper.GT.1.E-02) ) THEN - ! Defining mid-layer height (zl), interface height (zq), and layer depth (dz2). - ! This is already done above if the small-scale GWD is activated. - do k = kts,km - do i = its,im - zq(i,k+1) = PHII(i,k+1)*g_inv - dz2(i,k) = (PHII(i,k+1)-PHII(i,k))*g_inv - zl(i,k) = PHIL(i,k)*g_inv - enddo - enddo - ENDIF DO i=its,im IF ((xland(i)-1.5) .le. 0.) then diff --git a/physics/module_gfdl_cloud_microphys.F90 b/physics/module_gfdl_cloud_microphys.F90 index 2f6e5ec1a..01ab4655c 100644 --- a/physics/module_gfdl_cloud_microphys.F90 +++ b/physics/module_gfdl_cloud_microphys.F90 @@ -4729,7 +4729,7 @@ subroutine cloud_diagnosis (is, ie, ks, ke, den, delp, lsm, qmw, qmi, qmr, qms, real :: n0r = 8.0e6, n0s = 3.0e6, n0g = 4.0e6 real :: alphar = 0.8, alphas = 0.25, alphag = 0.5 real :: gammar = 17.837789, gammas = 8.2850630, gammag = 11.631769 - real :: qmin = 1.0e-12, beta = 1.22 + real :: qmin = 1.0e-12, beta = 1.22, qmin1 = 9.e-6 do k = ks, ke do i = is, ie @@ -4759,7 +4759,7 @@ subroutine cloud_diagnosis (is, ie, ks, ke, den, delp, lsm, qmw, qmi, qmr, qms, ! cloud ice (Heymsfield and Mcfarquhar, 1996) ! ----------------------------------------------------------------------- - if (qmi (i, k) .gt. qmin) then + if (qmi (i, k) .gt. qmin1) then qci (i, k) = dpg * qmi (i, k) * 1.0e3 rei_fac = log (1.0e3 * qmi (i, k) * den (i, k)) if (t (i, k) - tice .lt. - 50) then @@ -4785,7 +4785,7 @@ subroutine cloud_diagnosis (is, ie, ks, ke, den, delp, lsm, qmw, qmi, qmr, qms, ! cloud ice (Wyser, 1998) ! ----------------------------------------------------------------------- - if (qmi (i, k) .gt. qmin) then + if (qmi (i, k) .gt. qmin1) then qci (i, k) = dpg * qmi (i, k) * 1.0e3 bw = - 2. + 1.e-3 * log10 (den (i, k) * qmi (i, k) / rho_0) * max (0.0, tice - t (i, k)) ** 1.5 rei (i, k) = 377.4 + bw * (203.3 + bw * (37.91 + 2.3696 * bw)) @@ -4815,7 +4815,7 @@ subroutine cloud_diagnosis (is, ie, ks, ke, den, delp, lsm, qmw, qmi, qmr, qms, ! snow (Lin et al., 1983) ! ----------------------------------------------------------------------- - if (qms (i, k) .gt. qmin) then + if (qms (i, k) .gt. qmin1) then qcs (i, k) = dpg * qms (i, k) * 1.0e3 lambdas = exp (0.25 * log (pi * rhos * n0s / qms (i, k) / den (i, k))) res (i, k) = 0.5 * exp (log (gammas / 6) / alphas) / lambdas * 1.0e6 diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index 27552d9aa..b1ca6ba07 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -1302,7 +1302,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & kmax_qc = k qc_max = qc1d(k) elseif (qc1d(k) .lt. 0.0) then - write(*,*) 'WARNING, negative qc ', qc1d(k), & + write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qc ', qc1d(k), & ' at i,j,k=', i,j,k endif if (qr1d(k) .gt. qr_max) then @@ -1311,7 +1311,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & kmax_qr = k qr_max = qr1d(k) elseif (qr1d(k) .lt. 0.0) then - write(*,*) 'WARNING, negative qr ', qr1d(k), & + write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qr ', qr1d(k), & ' at i,j,k=', i,j,k endif if (nr1d(k) .gt. nr_max) then @@ -1320,7 +1320,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & kmax_nr = k nr_max = nr1d(k) elseif (nr1d(k) .lt. 0.0) then - write(*,*) 'WARNING, negative nr ', nr1d(k), & + write(*,'(a,e16.7,a,3i8)') 'WARNING, negative nr ', nr1d(k), & ' at i,j,k=', i,j,k endif if (qs1d(k) .gt. qs_max) then @@ -1329,7 +1329,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & kmax_qs = k qs_max = qs1d(k) elseif (qs1d(k) .lt. 0.0) then - write(*,*) 'WARNING, negative qs ', qs1d(k), & + write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qs ', qs1d(k), & ' at i,j,k=', i,j,k endif if (qi1d(k) .gt. qi_max) then @@ -1338,7 +1338,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & kmax_qi = k qi_max = qi1d(k) elseif (qi1d(k) .lt. 0.0) then - write(*,*) 'WARNING, negative qi ', qi1d(k), & + write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qi ', qi1d(k), & ' at i,j,k=', i,j,k endif if (qg1d(k) .gt. qg_max) then @@ -1347,7 +1347,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & kmax_qg = k qg_max = qg1d(k) elseif (qg1d(k) .lt. 0.0) then - write(*,*) 'WARNING, negative qg ', qg1d(k), & + write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qg ', qg1d(k), & ' at i,j,k=', i,j,k endif if (ni1d(k) .gt. ni_max) then @@ -1356,11 +1356,11 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & kmax_ni = k ni_max = ni1d(k) elseif (ni1d(k) .lt. 0.0) then - write(*,*) 'WARNING, negative ni ', ni1d(k), & + write(*,'(a,e16.7,a,3i8)') 'WARNING, negative ni ', ni1d(k), & ' at i,j,k=', i,j,k endif if (qv1d(k) .lt. 0.0) then - write(*,*) 'WARNING, negative qv ', qv1d(k), & + write(*,'(a,e16.7,a,3i8)') 'WARNING, negative qv ', qv1d(k), & ' at i,j,k=', i,j,k if (k.lt.kte-2 .and. k.gt.kts+1) then write(*,*) ' below and above are: ', qv(i,k-1,j), qv(i,k+1,j) diff --git a/physics/module_mp_thompson_make_number_concentrations.F90 b/physics/module_mp_thompson_make_number_concentrations.F90 index ef6779a67..b31753aa2 100644 --- a/physics/module_mp_thompson_make_number_concentrations.F90 +++ b/physics/module_mp_thompson_make_number_concentrations.F90 @@ -79,6 +79,11 @@ elemental real function make_IceNumber (Q_ice, temp) 161.503, 168.262, 175.248, 182.473, 189.952, 197.699, & 205.728, 214.055, 222.694, 231.661, 240.971, 250.639 /) + if (Q_ice == 0) then + make_IceNumber = 0 + return + end if + !+---+-----------------------------------------------------------------+ !..From the model 3D temperature field, subtract 179K for which !.. index value of retab as a start. Value of corr is for @@ -133,6 +138,11 @@ elemental real function make_DropletNumber (Q_cloud, qnwfa) real:: q_nwfa, x1, xDc integer:: nu_c + if (Q_cloud == 0) then + make_DropletNumber = 0 + return + end if + !+---+ q_nwfa = MAX(99.E6, MIN(qnwfa,5.E10)) @@ -160,6 +170,11 @@ elemental real function make_RainNumber (Q_rain, temp) !real, parameter:: PI = 3.1415926536 real, parameter:: am_r = PI*1000./6. + if (Q_rain == 0) then + make_RainNumber = 0 + return + end if + !+---+-----------------------------------------------------------------+ !.. Not thrilled with it, but set Y-intercept parameter to Marshal-Palmer value !.. that basically assumes melting snow becomes typical rain. However, for @@ -172,7 +187,7 @@ elemental real function make_RainNumber (Q_rain, temp) N0 = 8.E6 if (temp .le. 271.15) then - N0 = 8.E8 + N0 = 8.E8 elseif (temp .gt. 271.15 .and. temp.lt.273.15) then N0 = 8. * 10**(279.15-temp) endif diff --git a/physics/mp_thompson_post.F90 b/physics/mp_thompson_post.F90 index a21f668ec..feb031a3e 100644 --- a/physics/mp_thompson_post.F90 +++ b/physics/mp_thompson_post.F90 @@ -67,7 +67,7 @@ end subroutine mp_thompson_post_init !! #endif subroutine mp_thompson_post_run(ncol, nlev, tgrs_save, tgrs, prslk, dtp, & - mpicomm, mpirank, mpiroot, errmsg, errflg) + kdt, mpicomm, mpirank, mpiroot, errmsg, errflg) implicit none @@ -78,6 +78,7 @@ subroutine mp_thompson_post_run(ncol, nlev, tgrs_save, tgrs, prslk, dtp, & real(kind_phys), dimension(1:ncol,1:nlev), intent(inout) :: tgrs real(kind_phys), dimension(1:ncol,1:nlev), intent(in) :: prslk real(kind_phys), intent(in) :: dtp + integer, intent(in) :: kdt ! MPI information integer, intent(in ) :: mpicomm integer, intent(in ) :: mpirank @@ -115,8 +116,8 @@ subroutine mp_thompson_post_run(ncol, nlev, tgrs_save, tgrs, prslk, dtp, & if (tgrs_save(i,k) + mp_tend(i,k)*prslk(i,k) .ne. tgrs(i,k)) then #ifdef DEBUG - write(0,*) "mp_thompson_post_run mp_tend limiter: i, k, t_old, t_new, t_lim:", & - & i, k, tgrs_save(i,k), tgrs(i,k), tgrs_save(i,k) + mp_tend(i,k)*prslk(i,k) + write(0,'(a,3i6,3e16.7)') "mp_thompson_post_run mp_tend limiter: kdt, i, k, t_old, t_new, t_lim:", & + & kdt, i, k, tgrs_save(i,k), tgrs(i,k), tgrs_save(i,k) + mp_tend(i,k)*prslk(i,k) #endif events = events + 1 end if @@ -125,7 +126,8 @@ subroutine mp_thompson_post_run(ncol, nlev, tgrs_save, tgrs, prslk, dtp, & end do if (events > 0) then - write(0,'(a,i0,a,i0,a)') "mp_thompson_post_run: mp_tend_lim applied ", events, "/", nlev*ncol, " times" + write(0,'(a,i0,a,i0,a,i0)') "mp_thompson_post_run: mp_tend_lim applied ", events, "/", nlev*ncol, & + & " times at timestep ", kdt end if end subroutine mp_thompson_post_run diff --git a/physics/mp_thompson_post.meta b/physics/mp_thompson_post.meta index f1df2dd35..0f3cc6189 100644 --- a/physics/mp_thompson_post.meta +++ b/physics/mp_thompson_post.meta @@ -92,6 +92,14 @@ kind = kind_phys intent = in optional = F +[kdt] + standard_name = index_of_time_step + long_name = current forecast iteration + units = index + dimensions = () + type = integer + intent = in + optional = F [mpicomm] standard_name = mpi_comm long_name = MPI communicator diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index c3d061a9c..546cefca6 100644 --- a/physics/satmedmfvdifq.F +++ b/physics/satmedmfvdifq.F @@ -196,7 +196,7 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & & rlmn, rlmn1, rlmx, elmx, & ttend, utend, vtend, qtend, & zfac, zfmin, vk, spdk2, - & tkmin, xkzinv, xkgdx, + & tkmin, tkminx, xkzinv, xkgdx, & zlup, zldn, bsum, & tem, tem1, tem2, & ptem, ptem0, ptem1, ptem2 @@ -215,11 +215,11 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & parameter(prmin=0.25,prmax=4.0) parameter(pr0=1.0,prtke=1.0,prscu=0.67) parameter(f0=1.e-4,crbmin=0.15,crbmax=0.35) - parameter(tkmin=1.e-9,dspmax=10.0) + parameter(tkmin=1.e-9,tkminx=0.2,dspmax=10.0) parameter(qmin=1.e-8,qlmin=1.e-12,zfmin=1.e-8) parameter(aphi5=5.,aphi16=16.) parameter(elmfac=1.0,elefac=1.0,cql=100.) - parameter(dw2min=1.e-4,dkmax=1000.,xkgdx=25000.) + parameter(dw2min=1.e-4,dkmax=1000.,xkgdx=5000.) parameter(qlcr=3.5e-5,zstblmax=2500.,xkzinv=0.1) parameter(h1=0.33333333) parameter(ck0=0.4,ck1=0.15,ch0=0.4,ch1=0.15) @@ -326,20 +326,20 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & xkzo(i,k) = 0.0 xkzmo(i,k) = 0.0 if (k < kinver(i)) then -! vertical background diffusivity - ptem = prsi(i,k+1) * tx1(i) - tem1 = 1.0 - ptem - tem2 = tem1 * tem1 * 10.0 - tem2 = min(1.0, exp(-tem2)) - xkzo(i,k) = xkzm_hx(i) * tem2 -! +! minimum turbulent mixing length ptem = prsl(i,k) * tx1(i) tem1 = 1.0 - ptem tem2 = tem1 * tem1 * 2.5 tem2 = min(1.0, exp(-tem2)) rlmnz(i,k)= rlmn * tem2 rlmnz(i,k)= max(rlmnz(i,k), rlmn1) -! vertical background diffusivity for momentum +! vertical background diffusivity + ptem = prsi(i,k+1) * tx1(i) + tem1 = 1.0 - ptem + tem2 = tem1 * tem1 * 10.0 + tem2 = min(1.0, exp(-tem2)) + xkzo(i,k) = xkzm_hx(i) * tem2 +! vertical background diffusivity for momentum if (ptem >= xkzm_s) then xkzmo(i,k) = xkzm_mx(i) kx1(i) = k + 1 @@ -727,20 +727,20 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & ! ! background diffusivity decreasing with increasing surface layer stability ! - do i = 1, im - if(.not.sfcflg(i)) then - tem = (1. + 5. * rbsoil(i))**2. -! tem = (1. + 5. * zol(i))**2. - frik(i) = 0.1 + 0.9 / tem - endif - enddo -! - do k = 1,km1 - do i=1,im - xkzo(i,k) = frik(i) * xkzo(i,k) - xkzmo(i,k)= frik(i) * xkzmo(i,k) - enddo - enddo +! do i = 1, im +! if(.not.sfcflg(i)) then +! tem = (1. + 5. * rbsoil(i))**2. +!! tem = (1. + 5. * zol(i))**2. +! frik(i) = 0.1 + 0.9 / tem +! endif +! enddo +! +! do k = 1,km1 +! do i=1,im +! xkzo(i,k) = frik(i) * xkzo(i,k) +! xkzmo(i,k)= frik(i) * xkzmo(i,k) +! enddo +! enddo ! ! The background vertical diffusivities in the inversion layers are limited ! to be less than or equal to xkzminv @@ -920,13 +920,14 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & do i = 1, im if(k == 1) then tem = ckz(i,1) - tem1 = xkzmo(i,1) + tem1 = 0.5 * xkzmo(i,1) else tem = 0.5 * (ckz(i,k-1) + ckz(i,k)) tem1 = 0.5 * (xkzmo(i,k-1) + xkzmo(i,k)) endif ptem = tem1 / (tem * elm(i,k)) tkmnz(i,k) = ptem * ptem + tkmnz(i,k) = min(tkmnz(i,k), tkminx) tkmnz(i,k) = max(tkmnz(i,k), tkmin) enddo enddo