diff --git a/GFS_layer/GFS_physics_driver.F90 b/GFS_layer/GFS_physics_driver.F90 index 694d2191f..6e0c61ff8 100644 --- a/GFS_layer/GFS_physics_driver.F90 +++ b/GFS_layer/GFS_physics_driver.F90 @@ -20,6 +20,20 @@ module module_physics_driver GFS_control_type, GFS_grid_type, & GFS_tbd_type, GFS_cldprop_type, & GFS_radtend_type, GFS_diag_type + use edmf, only: edmf_run + use GFS_PBL_generic_pre, only: GFS_PBL_generic_pre_run + use GFS_PBL_generic_post, only: GFS_PBL_generic_post_run +! use sasas_deep, only: sasasdeep_run + use GFS_DCNV_generic_pre, only: GFS_DCNV_generic_pre_run + use GFS_DCNV_generic_post, only: GFS_DCNV_generic_post_run + use GFS_SCNV_generic_pre, only: GFS_SCNV_generic_pre_run + use GFS_SCNV_generic_post, only: GFS_SCNV_generic_post_run + use GFS_suite_interstitial_1, only: GFS_suite_interstitial_1_run + use GFS_suite_interstitial_2, only: GFS_suite_interstitial_2_run + use GFS_suite_interstitial_3, only: GFS_suite_interstitial_3_run + use GFS_suite_update_stateout, only: GFS_suite_update_stateout_run + use GFS_suite_interstitial_4, only: GFS_suite_interstitial_4_run + use GFS_suite_interstitial_5, only: GFS_suite_interstitial_5_run use GFS_zhaocarr_gscond, only: gscond_run use GFS_zhaocarr_precpd, only: precpd_run @@ -471,6 +485,9 @@ subroutine GFS_physics_driver & del, rhc, dtdt, dudt, dvdt, gwdcu, gwdcv, dtdtc, rainp, & ud_mf, dd_mf, dt_mf, prnum, dkt, sigmatot, sigmafrac + real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levs) :: & + initial_u, initial_v, initial_t, initial_qv + !--- GFDL modification for FV3 real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levs+1) ::& del_gz @@ -530,25 +547,26 @@ subroutine GFS_physics_driver & ! ! --- ... figure out number of extra tracers ! - tottracer = 0 ! no convective transport of tracers - if (Model%trans_trac .or. Model%cscnv) then - if (Model%ntcw > 0) then - if (Model%ntoz < Model%ntcw) then - trc_shft = Model%ntcw + Model%ncld - 1 - else - trc_shft = Model%ntoz - endif - elseif (Model%ntoz > 0) then - trc_shft = Model%ntoz - else - trc_shft = 1 - endif + ! tottracer = 0 ! no convective transport of tracers + ! if (Model%trans_trac .or. Model%cscnv) then + ! if (Model%ntcw > 0) then + ! if (Model%ntoz < Model%ntcw) then + ! trc_shft = Model%ntcw + Model%ncld - 1 + ! else + ! trc_shft = Model%ntoz + ! endif + ! elseif (Model%ntoz > 0) then + ! trc_shft = Model%ntoz + ! else + ! trc_shft = 1 + ! endif + ! + ! tracers = Model%ntrac - trc_shft + ! tottracer = tracers + ! if (Model%ntoz > 0) tottracer = tottracer + 1 ! ozone is added separately + ! endif + ! if (Model%ntke > 0) ntk = Model%ntke - trc_shft + 3 - tracers = Model%ntrac - trc_shft - tottracer = tracers - if (Model%ntoz > 0) tottracer = tottracer + 1 ! ozone is added separately - endif - if (Model%ntke > 0) ntk = Model%ntke - trc_shft + 3 ! if (lprnt) write(0,*)' trans_trac=',trans_trac,' tottracer=', & ! write(0,*)' trans_trac=',trans_trac,' tottracer=', & @@ -556,12 +574,14 @@ subroutine GFS_physics_driver & ! &, ntrac-ncld+2,' clstp=',clstp,' kdt=',kdt ! &,' ntk=',ntk,' lat=',lat - skip_macro = .false. + ! skip_macro = .false. + ! + ! allocate ( clw(ix,levs,tottracer+2) ) + ! if (Model%imfdeepcnv >= 0 .or. Model%imfshalcnv > 0) then + ! allocate (cnvc(ix,levs), cnvw(ix,levs)) + ! endif - allocate ( clw(ix,levs,tottracer+2) ) - if (Model%imfdeepcnv >= 0 .or. Model%imfshalcnv > 0) then - allocate (cnvc(ix,levs), cnvw(ix,levs)) - endif + call GFS_suite_interstitial_1_run (Model, Grid, tottracer, trc_shft, tracers, ntk, skip_macro, clw, cnvc, cnvw) ! ! --- set initial quantities for stochastic physics deltas if (Model%do_sppt) then @@ -606,19 +626,19 @@ subroutine GFS_physics_driver & call get_prs_fv3 (ix, levs, ntrac, Statein%phii, Statein%prsi, & Statein%tgrs, Statein%qgrs, del, del_gz) #endif + + call GFS_suite_interstitial_2_run (Model, Grid, Sfcprop, Statein, & + Diag, rhbbot, rhpbl, rhbtop, frain, islmsk, work1, work2, & + dudt, dvdt, dtdt, dtdtc, dqdt ) ! -!zhang: calrhc_run - rhbbot = Model%crtrh(1) - rhpbl = Model%crtrh(2) - rhbtop = Model%crtrh(3) ! ! --- ... frain=factor for centered difference scheme correction of rain amount. - frain = dtf / dtp + ! frain = dtf / dtp do i = 1, im sigmaf(i) = max( Sfcprop%vfrac(i),0.01 ) - islmsk(i) = nint(Sfcprop%slmsk(i)) + ! islmsk(i) = nint(Sfcprop%slmsk(i)) if (islmsk(i) == 2) then if (Model%isot == 1) then @@ -644,10 +664,10 @@ subroutine GFS_physics_driver & ! !GFDL work1(i) = (log(coslat(i) / (nlons(i)*latr)) - dxmin) * dxinv ! work1(i) = (log(Grid%dx(i)) - dxmin) * dxinv - work1(i) = (log(Grid%area(i)) - dxmin) * dxinv - work1(i) = max(0.0, min(1.0,work1(i))) - work2(i) = 1.0 - work1(i) - Diag%psurf(i) = Statein%pgr(i) + ! work1(i) = (log(Grid%area(i)) - dxmin) * dxinv + ! work1(i) = max(0.0, min(1.0,work1(i))) + ! work2(i) = 1.0 - work1(i) + ! Diag%psurf(i) = Statein%pgr(i) work3(i) = Statein%prsik(i,1) / Statein%prslk(i,1) !GFDL tem1 = con_rerth * (con_pi+con_pi)*coslat(i)/nlons(i) !GFDL tem2 = con_rerth * con_pi / latr @@ -681,11 +701,11 @@ subroutine GFS_physics_driver & smsoil(:,:) = Sfcprop%smc(:,:) stsoil(:,:) = Sfcprop%stc(:,:) slsoil(:,:) = Sfcprop%slc(:,:) !! clu: slc -> slsoil - dudt(:,:) = 0. - dvdt(:,:) = 0. - dtdt(:,:) = 0. - dtdtc(:,:) = 0. - dqdt(:,:,:) = 0. + ! dudt(:,:) = 0. + ! dvdt(:,:) = 0. + ! dtdt(:,:) = 0. + ! dtdtc(:,:) = 0. + ! dqdt(:,:,:) = 0. ! --- ... initialize dtdt with heating rate from dcyc2 @@ -761,51 +781,54 @@ subroutine GFS_physics_driver & gabsbdlw(:) = Radtend%semis(:) * adjsfcdlw(:) - if (Model%lssav) then ! --- ... accumulate/save output variables - -! --- ... sunshine duration time is defined as the length of time (in mdl output -! interval) that solar radiation falling on a plane perpendicular to the -! direction of the sun >= 120 w/m2 - - do i = 1, im - if ( xcosz(i) >= czmin ) then ! zenth angle > 89.994 deg - tem1 = adjsfcdsw(i) / xcosz(i) - if ( tem1 >= 120.0 ) then - Diag%suntim(i) = Diag%suntim(i) + dtf - endif - endif - enddo - -! --- ... sfc lw fluxes used by atmospheric model are saved for output - - if (Model%cplflx) then - do i = 1, im - if (flag_cice(i)) adjsfculw(i) = ulwsfc_cice(i) - enddo - endif - Diag%dlwsfc(:) = Diag%dlwsfc(:) + adjsfcdlw(:)*dtf - Diag%ulwsfc(:) = Diag%ulwsfc(:) + adjsfculw(:)*dtf - Diag%psmean(:) = Diag%psmean(:) + Statein%pgr(:)*dtf ! mean surface pressure - - if (Model%ldiag3d) then - if (Model%lsidea) then - Diag%dt3dt(:,:,1) = Diag%dt3dt(:,:,1) + Radtend%lwhd(:,:,1)*dtf - Diag%dt3dt(:,:,2) = Diag%dt3dt(:,:,2) + Radtend%lwhd(:,:,2)*dtf - Diag%dt3dt(:,:,3) = Diag%dt3dt(:,:,3) + Radtend%lwhd(:,:,3)*dtf - Diag%dt3dt(:,:,4) = Diag%dt3dt(:,:,4) + Radtend%lwhd(:,:,4)*dtf - Diag%dt3dt(:,:,5) = Diag%dt3dt(:,:,5) + Radtend%lwhd(:,:,5)*dtf - Diag%dt3dt(:,:,6) = Diag%dt3dt(:,:,6) + Radtend%lwhd(:,:,6)*dtf - else - do k = 1, levs - Diag%dt3dt(:,k,1) = Diag%dt3dt(:,k,1) + Radtend%htrlw(:,k)*dtf - Diag%dt3dt(:,k,2) = Diag%dt3dt(:,k,2) + Radtend%htrsw(:,k)*dtf*xmu(:) - enddo - endif - endif - endif ! end if_lssav_block +! if (Model%lssav) then ! --- ... accumulate/save output variables +! +! ! --- ... sunshine duration time is defined as the length of time (in mdl output +! ! interval) that solar radiation falling on a plane perpendicular to the +! ! direction of the sun >= 120 w/m2 +! +! do i = 1, im +! if ( xcosz(i) >= czmin ) then ! zenth angle > 89.994 deg +! tem1 = adjsfcdsw(i) / xcosz(i) +! if ( tem1 >= 120.0 ) then +! Diag%suntim(i) = Diag%suntim(i) + dtf +! endif +! endif +! enddo +! +! ! --- ... sfc lw fluxes used by atmospheric model are saved for output +! +! if (Model%cplflx) then +! do i = 1, im +! if (flag_cice(i)) adjsfculw(i) = ulwsfc_cice(i) +! enddo +! endif +! Diag%dlwsfc(:) = Diag%dlwsfc(:) + adjsfcdlw(:)*dtf +! Diag%ulwsfc(:) = Diag%ulwsfc(:) + adjsfculw(:)*dtf +! Diag%psmean(:) = Diag%psmean(:) + Statein%pgr(:)*dtf ! mean surface pressure +! +! if (Model%ldiag3d) then +! if (Model%lsidea) then +! Diag%dt3dt(:,:,1) = Diag%dt3dt(:,:,1) + Radtend%lwhd(:,:,1)*dtf +! Diag%dt3dt(:,:,2) = Diag%dt3dt(:,:,2) + Radtend%lwhd(:,:,2)*dtf +! Diag%dt3dt(:,:,3) = Diag%dt3dt(:,:,3) + Radtend%lwhd(:,:,3)*dtf +! Diag%dt3dt(:,:,4) = Diag%dt3dt(:,:,4) + Radtend%lwhd(:,:,4)*dtf +! Diag%dt3dt(:,:,5) = Diag%dt3dt(:,:,5) + Radtend%lwhd(:,:,5)*dtf +! Diag%dt3dt(:,:,6) = Diag%dt3dt(:,:,6) + Radtend%lwhd(:,:,6)*dtf +! else +! do k = 1, levs +! Diag%dt3dt(:,k,1) = Diag%dt3dt(:,k,1) + Radtend%htrlw(:,k)*dtf +! Diag%dt3dt(:,k,2) = Diag%dt3dt(:,k,2) + Radtend%htrsw(:,k)*dtf*xmu(:) +! enddo +! endif +! endif +! endif ! end if_lssav_block + call GFS_suite_interstitial_3_run (Model, Grid, Statein, Radtend, xcosz, & + adjsfcdsw, adjsfcdlw, adjsfculw, xmu, Diag, kcnv, hflx, evap) + call GFS_PBL_generic_pre_run (im, levs, kinver) - kcnv(:) = 0 - kinver(:) = levs + !kcnv(:) = 0 + !kinver(:) = levs invrsn(:) = .false. tx1(:) = 0.0 tx2(:) = 10.0 @@ -861,8 +884,8 @@ subroutine GFS_physics_driver & drain(:) = 0.0 ep1d(:) = 0.0 runof(:) = 0.0 - hflx(:) = 0.0 - evap(:) = 0.0 + !hflx(:) = 0.0 + !evap(:) = 0.0 evbs(:) = 0.0 evcw(:) = 0.0 trans(:) = 0.0 @@ -1060,10 +1083,10 @@ subroutine GFS_physics_driver & Diag%uswsfci(:) = adjsfcdsw(:) - adjsfcnsw(:) Diag%dswsfci(:) = adjsfcdsw(:) Diag%gfluxi(:) = gflx(:) - Diag%t1(:) = Statein%tgrs(:,1) - Diag%q1(:) = Statein%qgrs(:,1,1) - Diag%u1(:) = Statein%ugrs(:,1) - Diag%v1(:) = Statein%vgrs(:,1) + ! Diag%t1(:) = Statein%tgrs(:,1) + ! Diag%q1(:) = Statein%qgrs(:,1,1) + ! Diag%u1(:) = Statein%ugrs(:,1) + ! Diag%v1(:) = Statein%vgrs(:,1) ! --- ... update near surface fields @@ -1186,11 +1209,11 @@ subroutine GFS_physics_driver & Model%xkzm_m, Model%xkzm_h, Model%xkzm_s, lprnt, ipr, me) else if (Model%hybedmf) then - call moninedmf(ix, im, levs, nvdiff, Model%ntcw, dvdt, dudt, dtdt, dqdt,& + call edmf_run (ix, im, levs, nvdiff, Model%ntcw, dvdt, dudt, dtdt, dqdt,& Statein%ugrs, Statein%vgrs, Statein%tgrs, Statein%qgrs, & - Radtend%htrsw, Radtend%htrlw, xmu, Statein%prsik(1,1), & + Radtend%htrsw, Radtend%htrlw, xmu, Statein%prsik(:,1), & rb, Sfcprop%zorl, Diag%u10m, Diag%v10m, Sfcprop%ffmm, & - Sfcprop%ffhh, Sfcprop%tsfc, qss, hflx, evap, stress, & + Sfcprop%ffhh, Sfcprop%tsfc, hflx, evap, stress, & wind, kpbl, Statein%prsi, del, Statein%prsl, & Statein%prslk, Statein%phii, Statein%phil, dtp, & Model%dspheat, dusfc1, dvsfc1, dtsfc1, dqsfc1, Diag%hpbl,& @@ -1265,56 +1288,59 @@ subroutine GFS_physics_driver & Coupling%dtsfci_cpl(:) = dtsfc1(:) Coupling%dqsfci_cpl(:) = dqsfc1(:) endif -!-------------------------------------------------------lssav if loop ---------- - if (Model%lssav) then - Diag%dusfc (:) = Diag%dusfc(:) + dusfc1(:)*dtf - Diag%dvsfc (:) = Diag%dvsfc(:) + dvsfc1(:)*dtf - Diag%dtsfc (:) = Diag%dtsfc(:) + dtsfc1(:)*dtf - Diag%dqsfc (:) = Diag%dqsfc(:) + dqsfc1(:)*dtf - Diag%dusfci(:) = dusfc1(:) - Diag%dvsfci(:) = dvsfc1(:) - Diag%dtsfci(:) = dtsfc1(:) - Diag%dqsfci(:) = dqsfc1(:) -! if (lprnt) then -! write(0,*)' dusfc=',dusfc(ipr),' dusfc1=',dusfc1(ipr),' dtf=', -! & dtf,' kdt=',kdt,' lat=',lat -! endif - if (Model%ldiag3d) then - if (Model%lsidea) then - Diag%dt3dt(:,:,3) = Diag%dt3dt(:,:,3) + dtdt(:,:)*dtf - else - do k = 1, levs - do i = 1, im - tem = dtdt(i,k) - (Radtend%htrlw(i,k)+Radtend%htrsw(i,k)*xmu(i)) - Diag%dt3dt(i,k,3) = Diag%dt3dt(i,k,3) + tem*dtf - enddo - enddo - endif - Diag%du3dt(:,:,1) = Diag%du3dt(:,:,1) + dudt(:,:) * dtf - Diag%du3dt(:,:,2) = Diag%du3dt(:,:,2) - dudt(:,:) * dtf - Diag%dv3dt(:,:,1) = Diag%dv3dt(:,:,1) + dvdt(:,:) * dtf - Diag%dv3dt(:,:,2) = Diag%dv3dt(:,:,2) - dvdt(:,:) * dtf -! update dqdt_v to include moisture tendency due to vertical diffusion -! if (lgocart) then + call GFS_PBL_generic_post_run (Grid, Model, Radtend, dusfc1, dvsfc1, & + dtsfc1, dqsfc1, dudt, dvdt, dtdt, dqdt, xmu, Diag) +! !-------------------------------------------------------lssav if loop ---------- +! if (Model%lssav) then +! Diag%dusfc (:) = Diag%dusfc(:) + dusfc1(:)*dtf +! Diag%dvsfc (:) = Diag%dvsfc(:) + dvsfc1(:)*dtf +! Diag%dtsfc (:) = Diag%dtsfc(:) + dtsfc1(:)*dtf +! Diag%dqsfc (:) = Diag%dqsfc(:) + dqsfc1(:)*dtf +! Diag%dusfci(:) = dusfc1(:) +! Diag%dvsfci(:) = dvsfc1(:) +! Diag%dtsfci(:) = dtsfc1(:) +! Diag%dqsfci(:) = dqsfc1(:) +! ! if (lprnt) then +! ! write(0,*)' dusfc=',dusfc(ipr),' dusfc1=',dusfc1(ipr),' dtf=', +! ! & dtf,' kdt=',kdt,' lat=',lat +! ! endif +! +! if (Model%ldiag3d) then +! if (Model%lsidea) then +! Diag%dt3dt(:,:,3) = Diag%dt3dt(:,:,3) + dtdt(:,:)*dtf +! else +! do k = 1, levs +! do i = 1, im +! tem = dtdt(i,k) - (Radtend%htrlw(i,k)+Radtend%htrsw(i,k)*xmu(i)) +! Diag%dt3dt(i,k,3) = Diag%dt3dt(i,k,3) + tem*dtf +! enddo +! enddo +! endif +! Diag%du3dt(:,:,1) = Diag%du3dt(:,:,1) + dudt(:,:) * dtf +! Diag%du3dt(:,:,2) = Diag%du3dt(:,:,2) - dudt(:,:) * dtf +! Diag%dv3dt(:,:,1) = Diag%dv3dt(:,:,1) + dvdt(:,:) * dtf +! Diag%dv3dt(:,:,2) = Diag%dv3dt(:,:,2) - dvdt(:,:) * dtf +! ! update dqdt_v to include moisture tendency due to vertical diffusion +! ! if (lgocart) then +! ! do k = 1, levs +! ! do i = 1, im +! ! dqdt_v(i,k) = dqdt(i,k,1) * dtf +! ! enddo +! ! enddo +! ! endif ! do k = 1, levs ! do i = 1, im -! dqdt_v(i,k) = dqdt(i,k,1) * dtf +! tem = dqdt(i,k,1) * dtf +! Diag%dq3dt(i,k,1) = Diag%dq3dt(i,k,1) + tem ! enddo ! enddo +! if (Model%ntoz > 0) then +! Diag%dq3dt(:,:,5) = Diag%dq3dt(:,:,5) + dqdt(i,k,Model%ntoz) * dtf +! endif ! endif - do k = 1, levs - do i = 1, im - tem = dqdt(i,k,1) * dtf - Diag%dq3dt(i,k,1) = Diag%dq3dt(i,k,1) + tem - enddo - enddo - if (Model%ntoz > 0) then - Diag%dq3dt(:,:,5) = Diag%dq3dt(:,:,5) + dqdt(i,k,Model%ntoz) * dtf - endif - endif - - endif ! end if_lssav +! +! endif ! end if_lssav !-------------------------------------------------------lssav if loop ---------- ! ! Orographic gravity wave drag parameterization @@ -1400,10 +1426,12 @@ subroutine GFS_physics_driver & ! write(0,*)' dtdt=',(dtdt(ipr,ik),k=1,10) ! endif - Stateout%gt0(:,:) = Statein%tgrs(:,:) + dtdt(:,:) * dtp - Stateout%gu0(:,:) = Statein%ugrs(:,:) + dudt(:,:) * dtp - Stateout%gv0(:,:) = Statein%vgrs(:,:) + dvdt(:,:) * dtp - Stateout%gq0(:,:,:) = Statein%qgrs(:,:,:) + dqdt(:,:,:) * dtp + ! Stateout%gt0(:,:) = Statein%tgrs(:,:) + dtdt(:,:) * dtp + ! Stateout%gu0(:,:) = Statein%ugrs(:,:) + dudt(:,:) * dtp + ! Stateout%gv0(:,:) = Statein%vgrs(:,:) + dvdt(:,:) * dtp + ! Stateout%gq0(:,:,:) = Statein%qgrs(:,:,:) + dqdt(:,:,:) * dtp + + call GFS_suite_update_stateout_run (Statein, Model, Grid, dudt, dvdt, dtdt, dqdt, Stateout) ! if (lprnt) then ! write(7000,*)' ugrs=',ugrs(ipr,:) @@ -1490,17 +1518,19 @@ subroutine GFS_physics_driver & ! &,' lat=',lat,' kdt=',kdt,' me=',me ! if (lprnt) write(7000,*)' bef convection gv0=',gv0(ipr,:) - if (Model%ldiag3d) then - dtdt(:,:) = Stateout%gt0(:,:) - dudt(:,:) = Stateout%gu0(:,:) - dvdt(:,:) = Stateout%gv0(:,:) - elseif (Model%cnvgwd) then - dtdt(:,:) = Stateout%gt0(:,:) - endif ! end if_ldiag3d/cnvgwd + ! if (Model%ldiag3d) then + ! dtdt(:,:) = Stateout%gt0(:,:) + ! dudt(:,:) = Stateout%gu0(:,:) + ! dvdt(:,:) = Stateout%gv0(:,:) + ! elseif (Model%cnvgwd) then + ! dtdt(:,:) = Stateout%gt0(:,:) + ! endif ! end if_ldiag3d/cnvgwd + ! + ! if (Model%ldiag3d .or. Model%lgocart) then + ! dqdt(:,:,1) = Stateout%gq0(:,:,1) + ! endif ! end if_ldiag3d/lgocart - if (Model%ldiag3d .or. Model%lgocart) then - dqdt(:,:,1) = Stateout%gq0(:,:,1) - endif ! end if_ldiag3d/lgocart + call GFS_DCNV_generic_pre_run (Model, Stateout, Grid, initial_u, initial_v, initial_t, initial_qv) #ifdef GFS_HYDRO call get_phi(im, ix, levs, ntrac, Stateout%gt0, Stateout%gq0, & @@ -1517,13 +1547,14 @@ subroutine GFS_physics_driver & ! print *,' phii2=',phii(ipr,k=1,levs) ! print *,' phil2=',phil(ipr,:) ! endif - - clw(:,:,1) = 0.0 - clw(:,:,2) = -999.9 - if ((Model%imfdeepcnv >= 0) .or. (Model%imfshalcnv > 0)) then - cnvc(:,:) = 0.0 - cnvw(:,:) = 0.0 - endif + call GFS_suite_interstitial_4_run (Model, Grid, Statein, rhbbot, & + rhbtop, work1, work2, clw, cnvc, cnvw, ktop, kbot, rhc) + ! clw(:,:,1) = 0.0 + ! clw(:,:,2) = -999.9 + ! if ((Model%imfdeepcnv >= 0) .or. (Model%imfshalcnv > 0)) then + ! cnvc(:,:) = 0.0 + ! cnvw(:,:) = 0.0 + ! endif ! write(0,*)' before cnv clstp=',clstp,' kdt=',kdt,' lat=',lat @@ -1553,14 +1584,6 @@ subroutine GFS_physics_driver & ! -------------------------------------------- if (Model%ntcw > 0) then - do k=1,levs - do i=1,im -!zhang: gscond, precpd interstitial calrhc_run - tem = rhbbot - (rhbbot-rhbtop) * (1.0-Statein%prslk(i,k)) - tem = rhc_max * work1(i) + tem * work2(i) - rhc(i,k) = max(0.0, min(1.0,tem)) - enddo - enddo if (Model%ncld == 2) then clw(:,:,1) = Stateout%gq0(:,:,Model%ntiw) ! ice clw(:,:,2) = Stateout%gq0(:,:,Model%ntcw) ! water @@ -1858,45 +1881,48 @@ subroutine GFS_physics_driver & ! write(0,*)' aftcnvgq1=',(gq0(ipr,k,ntcw),k=1,levs) ! endif ! - do i = 1, im - Diag%rainc(:) = frain * rain1(:) - enddo +! do i = 1, im +! Diag%rainc(:) = frain * rain1(:) +! enddo +! ! +! if (Model%lssav) then +! Diag%cldwrk (:) = Diag%cldwrk (:) + cld1d(:) * dtf +! Diag%cnvprcp(:) = Diag%cnvprcp(:) + Diag%rainc(:) ! - if (Model%lssav) then - Diag%cldwrk (:) = Diag%cldwrk (:) + cld1d(:) * dtf - Diag%cnvprcp(:) = Diag%cnvprcp(:) + Diag%rainc(:) - - if (Model%ldiag3d) then - Diag%dt3dt(:,:,4) = Diag%dt3dt(:,:,4) + (Stateout%gt0(:,:)-dtdt(:,:)) * frain - Diag%dq3dt(:,:,2) = Diag%dq3dt(:,:,2) + (Stateout%gq0(:,:,1)-dqdt(:,:,1)) * frain - Diag%du3dt(:,:,3) = Diag%du3dt(:,:,3) + (Stateout%gu0(:,:)-dudt(:,:)) * frain - Diag%dv3dt(:,:,3) = Diag%dv3dt(:,:,3) + (Stateout%gv0(:,:)-dvdt(:,:)) * frain - - Diag%upd_mf(:,:) = Diag%upd_mf(:,:) + ud_mf(:,:) * (con_g*frain) - Diag%dwn_mf(:,:) = Diag%dwn_mf(:,:) + dd_mf(:,:) * (con_g*frain) - Diag%det_mf(:,:) = Diag%det_mf(:,:) + dt_mf(:,:) * (con_g*frain) - endif ! if (ldiag3d) +! if (Model%ldiag3d) then +! Diag%dt3dt(:,:,4) = Diag%dt3dt(:,:,4) + (Stateout%gt0(:,:)-dtdt(:,:)) * frain +! Diag%dq3dt(:,:,2) = Diag%dq3dt(:,:,2) + (Stateout%gq0(:,:,1)-dqdt(:,:,1)) * frain +! Diag%du3dt(:,:,3) = Diag%du3dt(:,:,3) + (Stateout%gu0(:,:)-dudt(:,:)) * frain +! Diag%dv3dt(:,:,3) = Diag%dv3dt(:,:,3) + (Stateout%gv0(:,:)-dvdt(:,:)) * frain +! +! Diag%upd_mf(:,:) = Diag%upd_mf(:,:) + ud_mf(:,:) * (con_g*frain) +! Diag%dwn_mf(:,:) = Diag%dwn_mf(:,:) + dd_mf(:,:) * (con_g*frain) +! Diag%det_mf(:,:) = Diag%det_mf(:,:) + dt_mf(:,:) * (con_g*frain) +! endif ! if (ldiag3d) +! +! endif ! end if_lssav - endif ! end if_lssav + call GFS_DCNV_generic_post_run (Grid, Model, Stateout, frain, rain1, cld1d, & + initial_u, initial_v, initial_t, initial_qv, ud_mf, dd_mf, dt_mf, cnvw, cnvc, Diag, Tbd) ! ! update dqdt_v to include moisture tendency due to deep convection if (Model%lgocart) then - Coupling%dqdti (:,:) = (Stateout%gq0(:,:,1) - dqdt(:,:,1)) * frain + Coupling%dqdti (:,:) = (Stateout%gq0(:,:,1) - initial_qv(:,:)) * frain Coupling%upd_mfi(:,:) = Coupling%upd_mfi(:,:) + ud_mf(:,:) * frain Coupling%dwn_mfi(:,:) = Coupling%dwn_mfi(:,:) + dd_mf(:,:) * frain Coupling%det_mfi(:,:) = Coupling%det_mfi(:,:) + dt_mf(:,:) * frain Coupling%cnvqci (:,:) = Coupling%cnvqci (:,:) + (clw(:,:,1)+clw(:,:,2))*frain endif ! if (lgocart) ! - if ((Model%npdf3d == 3) .and. (Model%num_p3d == 4)) then - num2 = Model%num_p3d + 2 - num3 = num2 + 1 - Tbd%phy_f3d(:,:,num2) = cnvw(:,:) - Tbd%phy_f3d(:,:,num3) = cnvc(:,:) - elseif ((Model%npdf3d == 0) .and. (Model%ncnvcld3d == 1)) then - num2 = Model%num_p3d + 1 - Tbd%phy_f3d(:,:,num2) = cnvw(:,:) - endif + ! if ((Model%npdf3d == 3) .and. (Model%num_p3d == 4)) then + ! num2 = Model%num_p3d + 2 + ! num3 = num2 + 1 + ! Tbd%phy_f3d(:,:,num2) = cnvw(:,:) + ! Tbd%phy_f3d(:,:,num3) = cnvc(:,:) + ! elseif ((Model%npdf3d == 0) .and. (Model%ncnvcld3d == 1)) then + ! num2 = Model%num_p3d + 1 + ! Tbd%phy_f3d(:,:,num2) = cnvw(:,:) + ! endif ! if (lprnt) write(7000,*)' bef cnvgwd gu0=',gu0(ipr,:) ! &,' lat=',lat,' kdt=',kdt,' me=',me @@ -1914,7 +1940,7 @@ subroutine GFS_physics_driver & do k = 1, levs do i = 1, im if (k >= kbot(i) .and. k <= ktop(i)) then - cumabs(i) = cumabs(i) + (Stateout%gt0(i,k)-dtdt(i,k)) * del(i,k) + cumabs(i) = cumabs(i) + (Stateout%gt0(i,k)-initial_t(i,k)) * del(i,k) work3(i) = work3(i) + del(i,k) endif enddo @@ -2071,12 +2097,14 @@ subroutine GFS_physics_driver & ! &,' lat=',lat,' kdt=',kdt,' me=',me !----------------Convective gravity wave drag parameterization over -------- - if (Model%ldiag3d) then - dtdt(:,:) = Stateout%gt0(:,:) - endif - if (Model%ldiag3d .or. Model%lgocart) then - dqdt(:,:,1) = Stateout%gq0(:,:,1) - endif + ! if (Model%ldiag3d) then + ! initial_t(:,:) = Stateout%gt0(:,:) + ! endif + ! if (Model%ldiag3d .or. Model%lgocart) then + ! initial_qv(:,:) = Stateout%gq0(:,:,1) + ! endif + + call GFS_SCNV_generic_pre_run (Model, Stateout, Grid, initial_t, initial_qv) ! write(0,*)' before do_shoc shal clstp=',clstp,' kdt=',kdt, ! & ' lat=',lat @@ -2157,27 +2185,29 @@ subroutine GFS_physics_driver & endif ! end if_imfshalcnv endif ! end if_shal_cnv - if (Model%lssav) then -! update dqdt_v to include moisture tendency due to shallow convection - if (Model%lgocart) then - do k = 1, levs - do i = 1, im - tem = (Stateout%gq0(i,k,1)-dqdt(i,k,1)) * frain - Coupling%dqdti(i,k) = Coupling%dqdti(i,k) + tem - enddo - enddo - endif - if (Model%ldiag3d) then - Diag%dt3dt(:,:,5) = Diag%dt3dt(:,:,5) + (Stateout%gt0(:,:)-dtdt(:,:)) * frain - Diag%dq3dt(:,:,3) = Diag%dq3dt(:,:,3) + (Stateout%gq0(:,:,1)-dqdt(:,:,1)) * frain - endif - endif ! end if_lssav -! - do k = 1, levs - do i = 1, im - if (clw(i,k,2) <= -999.0) clw(i,k,2) = 0.0 - enddo - enddo +! if (Model%lssav) then +! ! update dqdt_v to include moisture tendency due to shallow convection +! if (Model%lgocart) then +! do k = 1, levs +! do i = 1, im +! tem = (Stateout%gq0(i,k,1)-initial_qv(i,k)) * frain +! Coupling%dqdti(i,k) = Coupling%dqdti(i,k) + tem +! enddo +! enddo +! endif +! if (Model%ldiag3d) then +! Diag%dt3dt(:,:,5) = Diag%dt3dt(:,:,5) + (Stateout%gt0(:,:)-initial_t(:,:)) * frain +! Diag%dq3dt(:,:,3) = Diag%dq3dt(:,:,3) + (Stateout%gq0(:,:,1)-initial_qv(:,:)) * frain +! endif +! endif ! end if_lssav +! ! +! do k = 1, levs +! do i = 1, im +! if (clw(i,k,2) <= -999.0) clw(i,k,2) = 0.0 +! enddo +! enddo + + call GFS_SCNV_generic_post_run (Model, Stateout, Grid, initial_t, initial_qv, frain, Diag, clw) ! if (lprnt) then ! write(0,*)' prsl=',prsl(ipr,:) @@ -2317,26 +2347,17 @@ subroutine GFS_physics_driver & ! enddo ! endif if (Model%ldiag3d) then - Diag%dt3dt(:,:,4) = Diag%dt3dt(:,:,4) + (Stateout%gt0(:,:) -dtdt(:,:) ) * frain - Diag%dq3dt(:,:,2) = Diag%dq3dt(:,:,2) + (Stateout%gq0(:,:,1)-dqdt(:,:,1)) * frain + Diag%dt3dt(:,:,4) = Diag%dt3dt(:,:,4) + (Stateout%gt0(:,:) -initial_t(:,:) ) * frain + Diag%dq3dt(:,:,2) = Diag%dq3dt(:,:,2) + (Stateout%gq0(:,:,1)-initial_qv(:,:)) * frain endif endif endif ! moist convective adjustment over ! -!zhang -! if (Model%ldiag3d .or. Model%do_aw) then -! dtdt(:,:) = Stateout%gt0(:,:) -! dqdt(:,:,1) = Stateout%gq0(:,:,1) -! do n=Model%ntcw,Model%ntcw+Model%ncld-1 -! dqdt(:,:,n) = Stateout%gq0(:,:,n) -! enddo -! endif call GFS_MP_generic_pre_run (im, ix,levs,clw(:,:,1),clw(:,:,2), & Model%ldiag3d, Model%ntcw, Model%ncld, & Model%num_p3d, Stateout%gt0,Stateout%gq0(:,:,1), & dtdt,dqdt(:,:,1),dqdt(:,:,3) ) - ! dqdt_v : instaneous moisture tendency (kg/kg/sec) if (Model%lgocart) then Coupling%dqdti(:,:) = Coupling%dqdti(:,:) / dtf @@ -2541,8 +2562,8 @@ subroutine GFS_physics_driver & do k = 1,levs do i = 1,im tem1 = sigmafrac(i,k) - Stateout%gt0(i,k) = Stateout%gt0(i,k) - tem1 * (Stateout%gt0(i,k)-dtdt(i,k)) - tem2 = tem1 * (Stateout%gq0(i,k,1)-dqdt(i,k,1)) + Stateout%gt0(i,k) = Stateout%gt0(i,k) - tem1 * (Stateout%gt0(i,k)-initial_t(i,k)) + tem2 = tem1 * (Stateout%gq0(i,k,1)-initial_qv(i,k)) Stateout%gq0(i,k,1) = Stateout%gq0(i,k,1) - tem2 temrain1(i) = temrain1(i) - (Statein%prsi(i,k)-Statein%prsi(i,k+1)) & * tem2 * onebg @@ -2611,12 +2632,6 @@ subroutine GFS_physics_driver & ! if (Model%lssav) then ! Diag%totprcp(:) = Diag%totprcp(:) + Diag%rain(:) -! if (Model%ldiag3d) then -! Diag%dt3dt(:,:,6) = Diag%dt3dt(:,:,6) + (Stateout%gt0(:,:)-dtdt(:,:)) * frain -! Diag%dq3dt(:,:,4) = Diag%dq3dt(:,:,4) + (Stateout%gq0(:,:,1)-dqdt(:,:,1)) * frain -! endif -! endif - ! --- ... estimate t850 for rain-snow decision ! t850(:) = Stateout%gt0(:,1) @@ -2755,12 +2770,13 @@ subroutine GFS_physics_driver & enddo endif - deallocate (clw) + call GFS_suite_interstitial_5_run (clw, cnvc, cnvw) + !deallocate (clw) if (Model%do_shoc) then deallocate (qpl, qpi, ncpl, ncpi) endif - if (allocated(cnvc)) deallocate(cnvc) - if (allocated(cnvw)) deallocate(cnvw) + ! if (allocated(cnvc)) deallocate(cnvc) + ! if (allocated(cnvw)) deallocate(cnvw) ! deallocate (fscav, fswtr) ! diff --git a/makefile b/makefile index 16e10e7f3..c44cb823b 100644 --- a/makefile +++ b/makefile @@ -111,6 +111,7 @@ SRCS_f = \ ./physics/shalcv_fixdp.f \ ./physics/shalcv_opr.f \ ./physics/tracer_const_h.f \ + ./physics/tridi.f \ ./physics/tridi2t3.f SRCS_f90 = \ @@ -123,6 +124,10 @@ SRCS_f90 = \ ./physics/gcm_shoc.f90 \ ./physics/gcycle.f90 \ ./physics/get_prs_fv3.f90 \ + ./physics/GFS_DCNV_generic.f90 \ + ./physics/GFS_SCNV_generic.f90 \ + ./physics/GFS_PBL_generic.f90 \ + ./physics/GFS_suite_interstitial.f90 \ ./physics/h2ointerp.f90 \ ./physics/m_micro_driver.f90 \ ./physics/module_nst_model.f90 \ @@ -192,4 +197,3 @@ include ./depend ifneq (clean,$(findstring clean,$(MAKECMDGOALS))) -include depend endif - diff --git a/physics/GFS_DCNV_generic.f90 b/physics/GFS_DCNV_generic.f90 new file mode 100644 index 000000000..5891cbb78 --- /dev/null +++ b/physics/GFS_DCNV_generic.f90 @@ -0,0 +1,135 @@ +!> \file GFS_DCNV_generic.f90 +!! Contains code related to deep convective schemes to be used within the GFS physics suite. + + module GFS_DCNV_generic_pre + + contains + + subroutine GFS_DCNV_generic_pre_init () + end subroutine GFS_DCNV_generic_pre_init + + subroutine GFS_DCNV_generic_pre_finalize() + end subroutine GFS_DCNV_generic_pre_finalize + +!> \section arg_table_GFS_DCNV_generic_pre_run Argument Table +!! | local var name | longname | description | units | rank | type | kind | intent | optional | +!! |----------------|--------------------------------------------------------|-----------------------------------------------------------------------|---------------|------|-------------------------------|-----------|--------|----------| +!! | Model | FV3-GFS_Control_type | Fortran DDT containing FV3-GFS model control parameters | DDT | 0 | GFS_typedefs%GFS_control_type | | in | F | +!! | Stateout | FV3-GFS_Stateout_type | Fortran DDT containing FV3-GFS prognostic state to return to dycore | DDT | 0 | GFS_typedefs%GFS_stateout_type| | in | F | +!! | Grid | FV3-GFS_Grid_type | Fortran DDT containing FV3-GFS grid and interpolation related data | DDT | 0 | GFS_typedefs%GFS_grid_type | | in | F | +!! | initial_u | x_wind_initial | x-wind before entering a physics scheme | m s-1 | 2 | real | kind_phys | inout | F | +!! | initial_v | y_wind_initial | y-wind before entering a physics scheme | m s-1 | 2 | real | kind_phys | inout | F | +!! | initial_t | air_temperature_initial | air temperature before entering a physics scheme | K | 2 | real | kind_phys | inout | F | +!! | initial_qv | water_vapor_specific_humidity_initial | water vapor specific humidity before entering a physics scheme | kg kg-1 | 2 | real | kind_phys | inout | F | +!! + subroutine GFS_DCNV_generic_pre_run (Model, Stateout, Grid, initial_u, initial_v, initial_t, initial_qv) + use machine, only: kind_phys + use GFS_typedefs, only: GFS_control_type, GFS_stateout_type, GFS_grid_type + + type(GFS_control_type), intent(in) :: Model + type(GFS_stateout_type), intent(in) :: Stateout + type(GFS_grid_type), intent(in) :: Grid + real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levs), intent(inout) :: initial_u, initial_v, initial_t, initial_qv + + if (Model%ldiag3d) then + initial_t(:,:) = Stateout%gt0(:,:) + initial_u(:,:) = Stateout%gu0(:,:) + initial_v(:,:) = Stateout%gv0(:,:) + elseif (Model%cnvgwd) then + initial_t(:,:) = Stateout%gt0(:,:) + endif ! end if_ldiag3d/cnvgwd + + if (Model%ldiag3d .or. Model%lgocart) then + initial_qv(:,:) = Stateout%gq0(:,:,1) + endif ! end if_ldiag3d/lgocart + + end subroutine GFS_DCNV_generic_pre_run + + end module + + module GFS_DCNV_generic_post + + contains + + subroutine GFS_DCNV_generic_post_init () + end subroutine GFS_DCNV_generic_post_init + + subroutine GFS_DCNV_generic_post_finalize () + end subroutine GFS_DCNV_generic_post_finalize + +!> \section arg_table_GFS_DCNV_generic_post_run Argument Table +!! | local var name | longname | description | units | rank | type | kind | intent | optional | +!! |----------------|-----------------------------------------------------------|-----------------------------------------------------------------------|---------------|------|-------------------------------|-----------|--------|----------| +!! | Grid | FV3-GFS_Grid_type | Fortran DDT containing FV3-GFS grid and interpolation related data | DDT | 0 | GFS_typedefs%GFS_grid_type | | in | F | +!! | Model | FV3-GFS_Control_type | Fortran DDT containing FV3-GFS model control parameters | DDT | 0 | GFS_typedefs%GFS_control_type | | in | F | +!! | Stateout | FV3-GFS_Stateout_type | Fortran DDT containing FV3-GFS prognostic state to return to dycore | DDT | 0 | GFS_typedefs%GFS_stateout_type| | in | F | +!! | frain | dynamics_to_physics_timestep_ratio | ratio of dynamics timestep to physics timestep | none | 0 | real | kind_phys | in | F | +!! | rain1 | instantaneous_rainfall_amount | instantaneous rainfall amount | m | 1 | real | kind_phys | in | F | +!! | cld1d | cloud_work_function | cloud work function | m2 s-2 | 1 | real | kind_phys | in | F | +!! | initial_u | x_wind_initial | x-wind before entering a physics scheme | m s-1 | 2 | real | kind_phys | in | F | +!! | initial_v | y_wind_initial | y-wind before entering a physics scheme | m s-1 | 2 | real | kind_phys | in | F | +!! | initial_t | air_temperature_initial | air temperature before entering a physics scheme | K | 2 | real | kind_phys | in | F | +!! | initial_qv | water_vapor_specific_humidity_initial | water vapor specific humidity before entering a physics scheme | kg kg-1 | 2 | real | kind_phys | in | F | +!! | ud_mf | instantaneous_atmosphere_updraft_convective_mass_flux | (updraft mass flux) * delt | kg m-2 | 2 | real | kind_phys | in | F | +!! | dd_mf | instantaneous_atmosphere_downdraft_convective_mass_flux | (downdraft mass flux) * delt | kg m-2 | 2 | real | kind_phys | in | F | +!! | dt_mf | instantaneous_atmosphere_detrainment_convective_mass_flux | (detrainment mass flux) * delt | kg m-2 | 2 | real | kind_phys | in | F | +!! | cnvw | convective_cloud_water_specific_humidity | convective cloud water specific humidity | kg kg-1 | 2 | real | kind_phys | out | F | +!! | cnvc | convective_cloud_cover | convective cloud cover | frac | 2 | real | kind_phys | out | F | +!! | Diag | FV3-GFS_Diag_type | Fortran DDT containing FV3-GFS fields targeted for diagnostic output | DDT | 0 | GFS_typedefs%GFS_diag_type | | inout | F | +!! | Tbd | FV3-GFS_Tbd_type | Fortran DDT containing FV3-GFS miscellaneous data | DDT | 0 | GFS_typedefs%GFS_tbd_type | | inout | F | +!! + subroutine GFS_DCNV_generic_post_run (Grid, Model, Stateout, frain, rain1, cld1d, initial_u, initial_v, initial_t, initial_qv, & + ud_mf, dd_mf, dt_mf, cnvw, cnvc, Diag, Tbd) + + use machine, only: kind_phys + use GFS_typedefs, only: GFS_grid_type, GFS_control_type, GFS_stateout_type, GFS_diag_type, GFS_tbd_type + use physcons, only: con_g + + type(GFS_grid_type), intent(in) :: Grid + type(GFS_control_type), intent(in) :: Model + type(GFS_stateout_type), intent(in) :: Stateout + type(GFS_diag_type), intent(inout) :: Diag + type(GFS_tbd_type), intent(inout) :: Tbd + + real(kind=kind_phys), intent(in) :: frain + real(kind=kind_phys), dimension(size(Grid%xlon,1)), intent(in) :: rain1, cld1d + real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levs), intent(in) :: initial_u, initial_v, initial_t, initial_qv + real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levs), intent(in) :: ud_mf, dd_mf, dt_mf + real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levs), intent(in) :: cnvw, cnvc + + integer :: i, num2, num3 + + do i = 1, size(Grid%xlon,1) + Diag%rainc(:) = frain * rain1(:) + enddo + ! + if (Model%lssav) then + Diag%cldwrk (:) = Diag%cldwrk (:) + cld1d(:) * Model%dtf + Diag%cnvprcp(:) = Diag%cnvprcp(:) + Diag%rainc(:) + + if (Model%ldiag3d) then + Diag%dt3dt(:,:,4) = Diag%dt3dt(:,:,4) + (Stateout%gt0(:,:)-initial_t(:,:)) * frain + Diag%dq3dt(:,:,2) = Diag%dq3dt(:,:,2) + (Stateout%gq0(:,:,1)-initial_qv(:,:)) * frain + Diag%du3dt(:,:,3) = Diag%du3dt(:,:,3) + (Stateout%gu0(:,:)-initial_u(:,:)) * frain + Diag%dv3dt(:,:,3) = Diag%dv3dt(:,:,3) + (Stateout%gv0(:,:)-initial_v(:,:)) * frain + + Diag%upd_mf(:,:) = Diag%upd_mf(:,:) + ud_mf(:,:) * (con_g*frain) + Diag%dwn_mf(:,:) = Diag%dwn_mf(:,:) + dd_mf(:,:) * (con_g*frain) + Diag%det_mf(:,:) = Diag%det_mf(:,:) + dt_mf(:,:) * (con_g*frain) + endif ! if (ldiag3d) + + endif ! end if_lssav + + if ((Model%npdf3d == 3) .and. (Model%num_p3d == 4)) then + num2 = Model%num_p3d + 2 + num3 = num2 + 1 + Tbd%phy_f3d(:,:,num2) = cnvw(:,:) + Tbd%phy_f3d(:,:,num3) = cnvc(:,:) + elseif ((Model%npdf3d == 0) .and. (Model%ncnvcld3d == 1)) then + num2 = Model%num_p3d + 1 + Tbd%phy_f3d(:,:,num2) = cnvw(:,:) + endif + + end subroutine GFS_DCNV_generic_post_run + + end module diff --git a/physics/GFS_PBL_generic.f90 b/physics/GFS_PBL_generic.f90 new file mode 100644 index 000000000..de15fd0a4 --- /dev/null +++ b/physics/GFS_PBL_generic.f90 @@ -0,0 +1,126 @@ +!> \file GFS_PBL_generic.f90 +!! Contains code related to PBL schemes to be used within the GFS physics suite. + + module GFS_PBL_generic_pre + + contains + + subroutine GFS_PBL_generic_pre_init () + end subroutine GFS_PBL_generic_pre_init + + subroutine GFS_PBL_generic_pre_finalize() + end subroutine GFS_PBL_generic_pre_finalize + +!> \section arg_table_GFS_PBL_generic_pre_run Argument Table +!! | local var name | longname | description | units | rank | type | kind | intent | optional | +!! |----------------|--------------------------------------------------------|----------------------------------------------------|---------------|------|---------|-----------|--------|----------| +!! | im | horizontal_loop_extent | horizontal loop extent, start at 1 | index | 0 | integer | | in | F | +!! | levs | vertical_dimension | vertical layer dimension | index | 0 | integer | | in | F | +!! | kinver | index_of_highest_temperature_inversion | index of highest temperature inversion | index | 1 | integer | | in | F | +!! + subroutine GFS_PBL_generic_pre_run (im, levs, kinver) + + integer , intent(in) :: im, levs + integer, dimension(im), intent(inout) :: kinver + + kinver(:) = levs + + end subroutine GFS_PBL_generic_pre_run + + end module + + module GFS_PBL_generic_post + + contains + + subroutine GFS_PBL_generic_post_init () + end subroutine GFS_PBL_generic_post_init + + subroutine GFS_PBL_generic_post_finalize () + end subroutine GFS_PBL_generic_post_finalize + + + +!> \section arg_table_GFS_PBL_generic_post_run Argument Table +!! | local var name | longname | description | units | rank | type | kind | intent | optional | +!! |----------------|--------------------------------------------------------------|-----------------------------------------------------------------------|---------------|------|-------------------------------|-----------|--------|----------| +!! | Grid | FV3-GFS_Grid_type | Fortran DDT containing FV3-GFS grid and interpolation related data | DDT | 0 | GFS_typedefs%GFS_grid_type | | in | F | +!! | Model | FV3-GFS_Control_type | Fortran DDT containing FV3-GFS model control parameters | DDT | 0 | GFS_typedefs%GFS_control_type | | in | F | +!! | Radtend | FV3-GFS_Radtend_type | Fortran DDT containing FV3-GFS radiation tendencies needed in physics | DDT | 0 | GFS_typedefs%GFS_radtend_type | | in | F | +!! | dusfc1 | instantaneous_surface_x_momentum_flux | surface momentum flux in the x-direction valid for current call | Pa | 1 | real | kind_phys | in | F | +!! | dvsfc1 | instantaneous_surface_y_momentum_flux | surface momentum flux in the y-direction valid for current call | Pa | 1 | real | kind_phys | in | F | +!! | dtsfc1 | instantaneous_surface_upward_sensible_heat_flux | surface upward sensible heat flux valid for current call | W m-2 | 1 | real | kind_phys | in | F | +!! | dqsfc1 | instantaneous_surface_upward_latent_heat_flux | surface upward latent heat flux valid for current call | W m-2 | 1 | real | kind_phys | in | F | +!! | dudt | tendency_of_x_wind_due_to_model_physics | updated tendency of the x wind | m s-2 | 2 | real | kind_phys | in | F | +!! | dvdt | tendency_of_y_wind_due_to_model_physics | updated tendency of the y wind | m s-2 | 2 | real | kind_phys | in | F | +!! | dtdt | tendency_of_air_temperature_due_to_model_physics | updated tendency of the temperature | K s-1 | 2 | real | kind_phys | in | F | +!! | dqdt | tendency_of_tracers_due_to_model_physics | updated tendency of the tracers | kg kg-1 s-1 | 3 | real | kind_phys | in | F | +!! | xmu | zenith_angle_temporal_adjustment_factor_for_shortwave_fluxes | zenith angle temporal adjustment factor for shortwave fluxes | none | 2 | real | kind_phys | in | F | +!! | Diag | FV3-GFS_Diag_type | Fortran DDT containing FV3-GFS fields targeted for diagnostic output | DDT | 0 | GFS_typedefs%GFS_diag_type | | in | F | +!! + subroutine GFS_PBL_generic_post_run (Grid, Model, Radtend, dusfc1, dvsfc1, dtsfc1, dqsfc1, & + dudt, dvdt, dtdt, dqdt, xmu, Diag) + + use machine, only: kind_phys + use GFS_typedefs, only: GFS_diag_type, GFS_radtend_type, GFS_control_type, GFS_grid_type + + type(GFS_grid_type), intent(in) :: Grid + type(GFS_radtend_type), intent(in) :: Radtend + type(GFS_control_type), intent(in) :: Model + type(GFS_diag_type), intent(inout) :: Diag + + real(kind=kind_phys), dimension(size(Grid%xlon,1)), intent(in) :: dusfc1, dvsfc1, dtsfc1, dqsfc1, xmu + real(kind=kind_phys), dimension(size(Grid%xlon,1), Model%levs), intent(in) :: dudt, dvdt, dtdt + real(kind=kind_phys), dimension(size(Grid%xlon,1), Model%levs, Model%ntrac), intent(in) :: dqdt + + integer :: i, k + real(kind=kind_phys) :: tem + + if (Model%lssav) then + Diag%dusfc (:) = Diag%dusfc(:) + dusfc1(:)*Model%dtf + Diag%dvsfc (:) = Diag%dvsfc(:) + dvsfc1(:)*Model%dtf + Diag%dtsfc (:) = Diag%dtsfc(:) + dtsfc1(:)*Model%dtf + Diag%dqsfc (:) = Diag%dqsfc(:) + dqsfc1(:)*Model%dtf + Diag%dusfci(:) = dusfc1(:) + Diag%dvsfci(:) = dvsfc1(:) + Diag%dtsfci(:) = dtsfc1(:) + Diag%dqsfci(:) = dqsfc1(:) + ! if (lprnt) then + ! write(0,*)' dusfc=',dusfc(ipr),' dusfc1=',dusfc1(ipr),' dtf=', + ! & dtf,' kdt=',kdt,' lat=',lat + ! endif + + if (Model%ldiag3d) then + do k = 1, Model%levs + do i = 1, size(Grid%xlon,1) + tem = dtdt(i,k) - (Radtend%htrlw(i,k)+Radtend%htrsw(i,k)*xmu(i)) + Diag%dt3dt(i,k,3) = Diag%dt3dt(i,k,3) + tem*Model%dtf + enddo + enddo + Diag%du3dt(:,:,1) = Diag%du3dt(:,:,1) + dudt(:,:) * Model%dtf + Diag%du3dt(:,:,2) = Diag%du3dt(:,:,2) - dudt(:,:) * Model%dtf + Diag%dv3dt(:,:,1) = Diag%dv3dt(:,:,1) + dvdt(:,:) * Model%dtf + Diag%dv3dt(:,:,2) = Diag%dv3dt(:,:,2) - dvdt(:,:) * Model%dtf + ! update dqdt_v to include moisture tendency due to vertical diffusion + ! if (lgocart) then + ! do k = 1, levs + ! do i = 1, im + ! dqdt_v(i,k) = dqdt(i,k,1) * dtf + ! enddo + ! enddo + ! endif + do k = 1, Model%levs + do i = 1, size(Grid%xlon,1) + tem = dqdt(i,k,1) * Model%dtf + Diag%dq3dt(i,k,1) = Diag%dq3dt(i,k,1) + tem + enddo + enddo + if (Model%ntoz > 0) then + Diag%dq3dt(:,:,5) = Diag%dq3dt(:,:,5) + dqdt(i,k,Model%ntoz) * Model%dtf + endif + endif + + endif ! end if_lssav + end subroutine GFS_PBL_generic_post_run + + end module diff --git a/physics/GFS_SCNV_generic.f90 b/physics/GFS_SCNV_generic.f90 new file mode 100644 index 000000000..f504d892e --- /dev/null +++ b/physics/GFS_SCNV_generic.f90 @@ -0,0 +1,99 @@ +!> \file GFS_SCNV_generic.f90 +!! Contains code related to shallow convective schemes to be used within the GFS physics suite. + + module GFS_SCNV_generic_pre + + contains + + subroutine GFS_SCNV_generic_pre_init () + end subroutine GFS_SCNV_generic_pre_init + + subroutine GFS_SCNV_generic_pre_finalize() + end subroutine GFS_SCNV_generic_pre_finalize + +!> \section arg_table_GFS_SCNV_generic_pre_run Argument Table +!! | local var name | longname | description | units | rank | type | kind | intent | optional | +!! |----------------|--------------------------------------------------------|-----------------------------------------------------------------------|---------------|------|-------------------------------|-----------|--------|----------| +!! | Model | FV3-GFS_Control_type | Fortran DDT containing FV3-GFS model control parameters | DDT | 0 | GFS_typedefs%GFS_control_type | | in | F | +!! | Stateout | FV3-GFS_Stateout_type | Fortran DDT containing FV3-GFS prognostic state to return to dycore | DDT | 0 | GFS_typedefs%GFS_stateout_type| | in | F | +!! | Grid | FV3-GFS_Grid_type | Fortran DDT containing FV3-GFS grid and interpolation related data | DDT | 0 | GFS_typedefs%GFS_grid_type | | in | F | +!! | initial_t | air_temperature_initial | air temperature before entering a physics scheme | K | 2 | real | kind_phys | inout | F | +!! | initial_qv | water_vapor_specific_humidity_initial | water vapor specific humidity before entering a physics scheme | kg kg-1 | 2 | real | kind_phys | inout | F | +!! + subroutine GFS_SCNV_generic_pre_run (Model, Stateout, Grid, initial_t, initial_qv) + + use machine, only: kind_phys + use GFS_typedefs, only: GFS_control_type, GFS_stateout_type, GFS_grid_type + + type(GFS_control_type), intent(in) :: Model + type(GFS_stateout_type), intent(in) :: Stateout + type(GFS_grid_type), intent(in) :: Grid + real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levs), intent(inout) :: initial_t, initial_qv + + if (Model%ldiag3d) then + initial_t(:,:) = Stateout%gt0(:,:) + endif + if (Model%ldiag3d .or. Model%lgocart) then + initial_qv(:,:) = Stateout%gq0(:,:,1) + endif + + end subroutine GFS_SCNV_generic_pre_run + + end module + + module GFS_SCNV_generic_post + + contains + + subroutine GFS_SCNV_generic_post_init () + end subroutine GFS_SCNV_generic_post_init + + subroutine GFS_SCNV_generic_post_finalize () + end subroutine GFS_SCNV_generic_post_finalize + +!> \section arg_table_GFS_SCNV_generic_post_run Argument Table +!! | local var name | longname | description | units | rank | type | kind | intent | optional | +!! |----------------|-----------------------------------------------------------|-----------------------------------------------------------------------|---------------|------|-------------------------------|-----------|--------|----------| +!! | Model | FV3-GFS_Control_type | Fortran DDT containing FV3-GFS model control parameters | DDT | 0 | GFS_typedefs%GFS_control_type | | in | F | +!! | Stateout | FV3-GFS_Stateout_type | Fortran DDT containing FV3-GFS prognostic state to return to dycore | DDT | 0 | GFS_typedefs%GFS_stateout_type| | in | F | +!! | Grid | FV3-GFS_Grid_type | Fortran DDT containing FV3-GFS grid and interpolation related data | DDT | 0 | GFS_typedefs%GFS_grid_type | | in | F | +!! | initial_t | air_temperature_initial | air temperature before entering a physics scheme | K | 2 | real | kind_phys | in | F | +!! | initial_qv | water_vapor_specific_humidity_initial | water vapor specific humidity before entering a physics scheme | kg kg-1 | 2 | real | kind_phys | in | F | +!! | frain | dynamics_to_physics_timestep_ratio | ratio of dynamics timestep to physics timestep | none | 0 | real | kind_phys | in | F | +!! | Diag | FV3-GFS_Diag_type | Fortran DDT containing FV3-GFS fields targeted for diagnostic output | DDT | 0 | GFS_typedefs%GFS_diag_type | | inout | F | +!! | clw | convective_transportable_tracers | array to contain cloud water and other convective trans. tracers | kg kg-1 | 3 | real | kind_phys | inout | F | +!! + subroutine GFS_SCNV_generic_post_run (Model, Stateout, Grid, initial_t, initial_qv, frain, Diag, clw) + + use machine, only: kind_phys + use GFS_typedefs, only: GFS_grid_type, GFS_control_type, GFS_stateout_type, GFS_diag_type + + type(GFS_grid_type), intent(in) :: Grid + type(GFS_control_type), intent(in) :: Model + type(GFS_stateout_type), intent(in) :: Stateout + type(GFS_diag_type), intent(inout) :: Diag + + + real(kind=kind_phys), intent(in) :: frain + real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levs), intent(in) :: initial_t, initial_qv + + real(kind=kind_phys), intent(inout) :: clw(:,:,:) + + integer :: i, k + + if (Model%lssav) then + if (Model%ldiag3d) then + Diag%dt3dt(:,:,5) = Diag%dt3dt(:,:,5) + (Stateout%gt0(:,:)-initial_t(:,:)) * frain + Diag%dq3dt(:,:,3) = Diag%dq3dt(:,:,3) + (Stateout%gq0(:,:,1)-initial_qv(:,:)) * frain + endif + endif ! end if_lssav +! + do k = 1, Model%levs + do i = 1, size(Grid%xlon,1) + if (clw(i,k,2) <= -999.0) clw(i,k,2) = 0.0 + enddo + enddo + + end subroutine GFS_SCNV_generic_post_run + + end module diff --git a/physics/GFS_suite_interstitial.f90 b/physics/GFS_suite_interstitial.f90 new file mode 100644 index 000000000..d5c5ed67d --- /dev/null +++ b/physics/GFS_suite_interstitial.f90 @@ -0,0 +1,388 @@ +!> \file GFS_suite_interstitial.f90 +!! Contains code related to more than one scheme in the GFS physics suite. + + module GFS_suite_interstitial_1 + + contains + + subroutine GFS_suite_interstitial_1_init () + end subroutine GFS_suite_interstitial_1_init + + subroutine GFS_suite_interstitial_1_finalize() + end subroutine GFS_suite_interstitial_1_finalize + +!> \section arg_table_GFS_suite_interstitial_1_run Argument Table +!! | local var name | longname | description | units | rank | type | kind | intent | optional | +!! |----------------|--------------------------------------------------------|-----------------------------------------------------------------------|---------------|------|-------------------------------|-----------|--------|----------| +!! | Model | FV3-GFS_Control_type | Fortran DDT containing FV3-GFS model control parameters | DDT | 0 | GFS_typedefs%GFS_control_type | | in | F | +!! | Grid | FV3-GFS_Grid_type | Fortran DDT containing FV3-GFS grid and interpolation related data | DDT | 0 | GFS_typedefs%GFS_grid_type | | in | F | +!! | tottracer | number_of_total_tracers | total number of tracers | count | 0 | integer | | out | F | +!! | trc_shft | start_index_of_other_tracers | beginning index of the non-water tracer species | index | 0 | integer | | out | F | +!! | tracers | number_of_water_tracers | number of water-related tracers | index | 0 | integer | | out | F | +!! | ntk | index_of_TKE | index of TKE in the tracer array | index | 0 | integer | | out | F | +!! | skip_macro | flag_skip_macro | flag to skip cloud macrophysics in Morrison scheme | flag | 1 | logical | | out | F | +!! | clw | convective_transportable_tracers | array to contain cloud water and other convective trans. tracers | kg kg-1 | 3 | real | kind_phys | out | F | +!! | cnvc | convective_cloud_cover | convective cloud cover | frac | 2 | real | kind_phys | out | F | +!! | cnvw | convective_cloud_water_specific_humidity | convective cloud water specific humidity | kg kg-1 | 2 | real | kind_phys | out | F | +!! + subroutine GFS_suite_interstitial_1_run (Model, Grid, tottracer, trc_shft, tracers, ntk, skip_macro, clw, cnvc, cnvw) + + use machine, only: kind_phys + use GFS_typedefs, only: GFS_control_type, GFS_grid_type + + type(GFS_control_type), intent(in) :: Model + type(GFS_grid_type), intent(in) :: Grid + integer, intent(out) :: tottracer, trc_shft, tracers, ntk + logical, dimension(size(Grid%xlon,1)), intent(out) :: skip_macro + real(kind=kind_phys), allocatable, intent(out) :: clw(:,:,:), cnvc(:,:), cnvw(:,:) + + tottracer = 0 ! no convective transport of tracers + if (Model%trans_trac .or. Model%cscnv) then + if (Model%ntcw > 0) then + if (Model%ntoz < Model%ntcw) then + trc_shft = Model%ntcw + Model%ncld - 1 + else + trc_shft = Model%ntoz + endif + elseif (Model%ntoz > 0) then + trc_shft = Model%ntoz + else + trc_shft = 1 + endif + + tracers = Model%ntrac - trc_shft + tottracer = tracers + if (Model%ntoz > 0) tottracer = tottracer + 1 ! ozone is added separately + endif + if (Model%ntke > 0) ntk = Model%ntke - trc_shft + 3 + + skip_macro = .false. + + allocate ( clw(size(Grid%xlon,1),Model%levs,tottracer+2) ) + if (Model%imfdeepcnv >= 0 .or. Model%imfshalcnv > 0) then + allocate (cnvc(size(Grid%xlon,1),Model%levs), cnvw(size(Grid%xlon,1),Model%levs)) + endif + + end subroutine GFS_suite_interstitial_1_run + + end module + + module GFS_suite_interstitial_2 + + contains + + subroutine GFS_suite_interstitial_2_init () + end subroutine GFS_suite_interstitial_2_init + + subroutine GFS_suite_interstitial_2_finalize() + end subroutine GFS_suite_interstitial_2_finalize + +!> \section arg_table_GFS_suite_interstitial_2_run Argument Table +!! | local var name | longname | description | units | rank | type | kind | intent | optional | +!! |----------------|-------------------------------------------------------------------------|-------------------------------------------------------------------------|---------------|------|-------------------------------|-----------|--------|----------| +!! | Model | FV3-GFS_Control_type | Fortran DDT containing FV3-GFS model control parameters | DDT | 0 | GFS_typedefs%GFS_control_type | | in | F | +!! | Grid | FV3-GFS_Grid_type | Fortran DDT containing FV3-GFS grid and interpolation related data | DDT | 0 | GFS_typedefs%GFS_grid_type | | in | F | +!! | Sfcprop | FV3-GFS_Sfcprop_type | Fortran DDT containing FV3-GFS surface fields | DDT | 0 | GFS_typedefs%GFS_sfcprop_type | | in | F | +!! | Statein | FV3-GFS_Statein_type | Fortran DDT containing FV3-GFS prognostic state data in from dycore | DDT | 0 | GFS_typedefs%GFS_statein_type | | in | F | +!! | Diag | FV3-GFS_Diag_type | Fortran DDT containing FV3-GFS fields targeted for diagnostic output | DDT | 0 | GFS_typedefs%GFS_diag_type | | inout | F | +!! | rhbbot | critical_relative_humidity_at_surface | critical relative humidity at the surface | frac | 0 | real | kind_phys | out | F | +!! | rhpbl | critical_relative_humidity_at_PBL_top | critical relative humidity at the PBL top | frac | 0 | real | kind_phys | out | F | +!! | rhbtop | critical_relative_humidity_at_top_of_atmosphere | critical relative humidity at the top of atmosphere | frac | 0 | real | kind_phys | out | F | +!! | frain | dynamics_to_physics_timestep_ratio | ratio of dynamics timestep to physics timestep | none | 0 | real | kind_phys | out | F | +!! | islmsk | sea_land_ice_mask | landmask: sea/land/ice=0/1/2 | flag | 1 | integer | | out | F | +!! | work1 | grid_related_coefficient | grid size related coefficient used in scale-sensitive schemes | none | 1 | real | kind_phys | out | F | +!! | work2 | grid_related_coefficient_complement | complement to work1 | none | 1 | real | kind_phys | out | F | +!! | dudt | tendency_of_x_wind_due_to_model_physics | updated tendency of the x wind | m s-2 | 2 | real | kind_phys | out | F | +!! | dvdt | tendency_of_y_wind_due_to_model_physics | updated tendency of the y wind | m s-2 | 2 | real | kind_phys | out | F | +!! | dtdt | tendency_of_air_temperature_due_to_model_physics | updated tendency of the temperature | K s-1 | 2 | real | kind_phys | out | F | +!! | dtdtc | tendency_of_air_temperature_due_to_radiative_heating_assuming_clear_sky | clear sky radiative (shortwave + longwave) heating rate at current time | K s-1 | 2 | real | kind_phys | out | F | +!! | dqdt | tendency_of_tracers_due_to_model_physics | updated tendency of the tracers | kg kg-1 s-1 | 3 | real | kind_phys | out | F | +!! + subroutine GFS_suite_interstitial_2_run (Model, Grid, Sfcprop, Statein, Diag, rhbbot, rhpbl, rhbtop, frain, islmsk, work1, work2, dudt, dvdt, dtdt, dtdtc, dqdt) + + use machine, only: kind_phys + use physcons, only: dxmin, dxinv + use GFS_typedefs, only: GFS_control_type, GFS_grid_type, GFS_sfcprop_type, GFS_statein_type, GFS_diag_type + + type(GFS_control_type), intent(in) :: Model + type(GFS_grid_type), intent(in) :: Grid + type(GFS_sfcprop_type), intent(in) :: Sfcprop + type(GFS_statein_type), intent(in) :: Statein + type(GFS_diag_type), intent(inout) :: Diag + + real(kind=kind_phys), intent(out) :: rhbbot, rhpbl, rhbtop, frain + integer, dimension(size(Grid%xlon,1)), intent(out) :: islmsk + real(kind=kind_phys), dimension(size(Grid%xlon,1)), intent(out) :: work1, work2 + real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levs), intent(out) :: dudt, dvdt, dtdt, dtdtc + real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levs,Model%ntrac), intent(out) :: dqdt + + integer :: i + + rhbbot = Model%crtrh(1) + rhpbl = Model%crtrh(2) + rhbtop = Model%crtrh(3) + + frain = Model%dtf / Model%dtp + + do i = 1, size(Grid%xlon,1) + islmsk(i) = nint(Sfcprop%slmsk(i)) + work1(i) = (log(Grid%area(i)) - dxmin) * dxinv + work1(i) = max(0.0, min(1.0,work1(i))) + work2(i) = 1.0 - work1(i) + Diag%psurf(i) = Statein%pgr(i) + end do + + dudt(:,:) = 0. + dvdt(:,:) = 0. + dtdt(:,:) = 0. + dtdtc(:,:) = 0. + dqdt(:,:,:) = 0. + + end subroutine GFS_suite_interstitial_2_run + + end module + + module GFS_suite_interstitial_3 + + contains + + subroutine GFS_suite_interstitial_3_init () + end subroutine GFS_suite_interstitial_3_init + + subroutine GFS_suite_interstitial_3_finalize() + end subroutine GFS_suite_interstitial_3_finalize + +!> \section arg_table_GFS_suite_interstitial_3_run Argument Table +!! | local var name | longname | description | units | rank | type | kind | intent | optional | +!! |----------------|--------------------------------------------------------------|-----------------------------------------------------------------------|---------------|------|-------------------------------|-----------|--------|----------| +!! | Model | FV3-GFS_Control_type | Fortran DDT containing FV3-GFS model control parameters | DDT | 0 | GFS_typedefs%GFS_control_type | | in | F | +!! | Grid | FV3-GFS_Grid_type | Fortran DDT containing FV3-GFS grid and interpolation related data | DDT | 0 | GFS_typedefs%GFS_grid_type | | in | F | +!! | Statein | FV3-GFS_Statein_type | Fortran DDT containing FV3-GFS prognostic state data in from dycore | DDT | 0 | GFS_typedefs%GFS_statein_type | | in | F | +!! | Radtend | FV3-GFS_Radtend_type | Fortran DDT containing FV3-GFS radiation tendencies needed in physics | DDT | 0 | GFS_typedefs%GFS_radtend_type | | in | F | +!! | xcosz | instantaneous_cosine_of_zenith_angle | cosine of zenith angle at current time | none | 1 | real | kind_phys | in | F | +!! | adjsfcdsw | surface_downwelling_shortwave_flux | surface downwelling shortwave flux at current time | W m-2 | 1 | real | kind_phys | in | F | +!! | adjsfcdlw | surface_downwelling_longwave_flux | surface downwelling longwave flux at current time | W m-2 | 1 | real | kind_phys | in | F | +!! | adjsfculw | surface_upwelling_longwave_flux | surface upwelling longwave flux at current time | W m-2 | 1 | real | kind_phys | in | F | +!! | xmu | zenith_angle_temporal_adjustment_factor_for_shortwave_fluxes | zenith angle temporal adjustment factor for shortwave fluxes | none | 1 | real | kind_phys | in | F | +!! | Diag | FV3-GFS_Diag_type | Fortran DDT containing FV3-GFS fields targeted for diagnostic output | DDT | 0 | GFS_typedefs%GFS_diag_type | | inout | F | +!! | kcnv | flag_deep_convection | flag indicating whether convection occurs in column (0 or 1) | index | 1 | integer | | out | F | +!! | heat | kinematic_surface_upward_sensible_heat_flux | kinematic surface upward sensible heat flux | K m s-1 | 1 | real | kind_phys | out | F | +!! | evap | kinematic_surface_upward_latent_heat_flux | kinematic surface upward latent heat flux | kg kg-1 m s-1 | 1 | real | kind_phys | out | F | +!! + subroutine GFS_suite_interstitial_3_run (Model, Grid, Statein, Radtend, xcosz, adjsfcdsw, adjsfcdlw, adjsfculw, xmu, Diag, kcnv, hflx, evap) + + use machine, only: kind_phys + use GFS_typedefs, only: GFS_control_type, GFS_grid_type, GFS_statein_type, GFS_radtend_type, GFS_diag_type + + type(GFS_control_type), intent(in) :: Model + type(GFS_grid_type), intent(in) :: Grid + type(GFS_statein_type), intent(in) :: Statein + type(GFS_radtend_type), intent(in) :: Radtend + type(GFS_diag_type), intent(inout) :: Diag + + integer, dimension(size(Grid%xlon,1)), intent(out) :: kcnv + real(kind=kind_phys), dimension(size(Grid%xlon,1)), intent(in) :: xcosz, adjsfcdsw, adjsfcdlw, adjsfculw, xmu + real(kind=kind_phys), dimension(size(Grid%xlon,1)), intent(out) :: hflx, evap + + real(kind=kind_phys), parameter :: czmin = 0.0001 ! cos(89.994) + + integer :: i + + real(kind=kind_phys) :: tem1 + + if (Model%lssav) then ! --- ... accumulate/save output variables + + ! --- ... sunshine duration time is defined as the length of time (in mdl output + ! interval) that solar radiation falling on a plane perpendicular to the + ! direction of the sun >= 120 w/m2 + + do i = 1, size(Grid%xlon,1) + if ( xcosz(i) >= czmin ) then ! zenth angle > 89.994 deg + tem1 = adjsfcdsw(i) / xcosz(i) + if ( tem1 >= 120.0 ) then + Diag%suntim(i) = Diag%suntim(i) + Model%dtf + endif + endif + enddo + + ! --- ... sfc lw fluxes used by atmospheric model are saved for output + + Diag%dlwsfc(:) = Diag%dlwsfc(:) + adjsfcdlw(:)*Model%dtf + Diag%ulwsfc(:) = Diag%ulwsfc(:) + adjsfculw(:)*Model%dtf + Diag%psmean(:) = Diag%psmean(:) + Statein%pgr(:)*Model%dtf ! mean surface pressure + + if (Model%ldiag3d) then + if (Model%lsidea) then + Diag%dt3dt(:,:,1) = Diag%dt3dt(:,:,1) + Radtend%lwhd(:,:,1)*Model%dtf + Diag%dt3dt(:,:,2) = Diag%dt3dt(:,:,2) + Radtend%lwhd(:,:,2)*Model%dtf + Diag%dt3dt(:,:,3) = Diag%dt3dt(:,:,3) + Radtend%lwhd(:,:,3)*Model%dtf + Diag%dt3dt(:,:,4) = Diag%dt3dt(:,:,4) + Radtend%lwhd(:,:,4)*Model%dtf + Diag%dt3dt(:,:,5) = Diag%dt3dt(:,:,5) + Radtend%lwhd(:,:,5)*Model%dtf + Diag%dt3dt(:,:,6) = Diag%dt3dt(:,:,6) + Radtend%lwhd(:,:,6)*Model%dtf + else + do k = 1, Model%levs + Diag%dt3dt(:,k,1) = Diag%dt3dt(:,k,1) + Radtend%htrlw(:,k)*Model%dtf + Diag%dt3dt(:,k,2) = Diag%dt3dt(:,k,2) + Radtend%htrsw(:,k)*Model%dtf*xmu(:) + enddo + endif + endif + endif ! end if_lssav_block + + kcnv(:) = 0 + + hflx(:) = 0.0 + evap(:) = 0.0 + + Diag%t1(:) = Statein%tgrs(:,1) + Diag%q1(:) = Statein%qgrs(:,1,1) + Diag%u1(:) = Statein%ugrs(:,1) + Diag%v1(:) = Statein%vgrs(:,1) + + end subroutine GFS_suite_interstitial_3_run + + end module + + module GFS_suite_update_stateout + + contains + + subroutine GFS_suite_update_stateout_init () + end subroutine GFS_suite_update_stateout_init + + subroutine GFS_suite_update_stateout_finalize() + end subroutine GFS_suite_update_stateout_finalize + +!> \section arg_table_GFS_suite_update_stateout_run Argument Table +!! | local var name | longname | description | units | rank | type | kind | intent | optional | +!! |----------------|--------------------------------------------------------|-----------------------------------------------------------------------|---------------|------|-------------------------------|-----------|--------|----------| +!! | Statein | FV3-GFS_Statein_type | Fortran DDT containing FV3-GFS prognostic state data in from dycore | DDT | 0 | GFS_typedefs%GFS_statein_type | | in | F | +!! | Model | FV3-GFS_Control_type | Fortran DDT containing FV3-GFS model control parameters | DDT | 0 | GFS_typedefs%GFS_control_type | | in | F | +!! | Grid | FV3-GFS_Grid_type | Fortran DDT containing FV3-GFS grid and interpolation related data | DDT | 0 | GFS_typedefs%GFS_grid_type | | in | F | +!! | dudt | tendency_of_x_wind_due_to_model_physics | updated tendency of the x wind | m s-2 | 2 | real | kind_phys | in | F | +!! | dvdt | tendency_of_y_wind_due_to_model_physics | updated tendency of the y wind | m s-2 | 2 | real | kind_phys | in | F | +!! | dtdt | tendency_of_air_temperature_due_to_model_physics | updated tendency of the temperature | K s-1 | 2 | real | kind_phys | in | F | +!! | dqdt | tendency_of_tracers_due_to_model_physics | updated tendency of the tracers | kg kg-1 s-1 | 3 | real | kind_phys | in | F | +!! | Stateout | FV3-GFS_Stateout_type | Fortran DDT containing FV3-GFS prognostic state to return to dycore | DDT | 0 | GFS_typedefs%GFS_stateout_type| | inout | F | +!! + subroutine GFS_suite_update_stateout_run (Statein, Model, Grid, dudt, dvdt, dtdt, dqdt, Stateout) + + use machine, only: kind_phys + use GFS_typedefs, only: GFS_control_type, GFS_statein_type, GFS_grid_type, GFS_stateout_type + + type(GFS_control_type), intent(in) :: Model + type(GFS_statein_type), intent(in) :: Statein + type(GFS_grid_type), intent(in) :: Grid + type(GFS_stateout_type), intent(inout) :: Stateout + + real(kind=kind_phys), dimension(size(Grid%xlon,1), Model%levs), intent(in) :: dudt, dvdt, dtdt + real(kind=kind_phys), dimension(size(Grid%xlon,1), Model%levs, Model%ntrac), intent(in) :: dqdt + + Stateout%gt0(:,:) = Statein%tgrs(:,:) + dtdt(:,:) * Model%dtp + Stateout%gu0(:,:) = Statein%ugrs(:,:) + dudt(:,:) * Model%dtp + Stateout%gv0(:,:) = Statein%vgrs(:,:) + dvdt(:,:) * Model%dtp + Stateout%gq0(:,:,:) = Statein%qgrs(:,:,:) + dqdt(:,:,:) * Model%dtp + + end subroutine GFS_suite_update_stateout_run + +end module + +module GFS_suite_interstitial_4 + +contains + +subroutine GFS_suite_interstitial_4_init () +end subroutine GFS_suite_interstitial_4_init + +subroutine GFS_suite_interstitial_4_finalize() +end subroutine GFS_suite_interstitial_4_finalize + +!> \section arg_table_GFS_suite_interstitial_4_run Argument Table +!! | local var name | longname | description | units | rank | type | kind | intent | optional | +!! |----------------|--------------------------------------------------------|-----------------------------------------------------------------------|---------------|------|-------------------------------|-----------|--------|----------| +!! | Model | FV3-GFS_Control_type | Fortran DDT containing FV3-GFS model control parameters | DDT | 0 | GFS_typedefs%GFS_control_type | | in | F | +!! | Grid | FV3-GFS_Grid_type | Fortran DDT containing FV3-GFS grid and interpolation related data | DDT | 0 | GFS_typedefs%GFS_grid_type | | in | F | +!! | Statein | FV3-GFS_Statein_type | Fortran DDT containing FV3-GFS prognostic state data in from dycore | DDT | 0 | GFS_typedefs%GFS_statein_type | | in | F | +!! | rhbbot | critical_relative_humidity_at_surface | critical relative humidity at the surface | frac | 0 | real | kind_phys | in | F | +!! | rhbtop | critical_relative_humidity_at_top_of_atmosphere | critical relative humidity at the top of atmosphere | frac | 0 | real | kind_phys | in | F | +!! | work1 | grid_related_coefficient | grid size related coefficient used in scale-sensitive schemes | none | 1 | real | kind_phys | in | F | +!! | work2 | grid_related_coefficient_complement | complement to work1 | none | 1 | real | kind_phys | in | F | +!! | clw | convective_transportable_tracers | array to contain cloud water and other convective trans. tracers | kg kg-1 | 3 | real | kind_phys | inout | F | +!! | cnvc | convective_cloud_cover | convective cloud cover | frac | 2 | real | kind_phys | inout | F | +!! | cnvw | convective_cloud_water_specific_humidity | convective cloud water specific humidity | kg kg-1 | 2 | real | kind_phys | inout | F | +!! | ktop | vertical_index_at_cloud_top | vertical index at cloud top | index | 1 | integer | | inout | F | +!! | kbot | vertical_index_at_cloud_base | vertical index at cloud base | index | 1 | integer | | inout | F | +!! | rhc | critical_relative_humidity | critical relative humidity | frac | 2 | real | kind_phys | out | F | +!! +subroutine GFS_suite_interstitial_4_run (Model, Grid, Statein, rhbbot, rhbtop, work1, work2, clw, cnvc, cnvw, ktop, kbot, rhc) + + use machine, only: kind_phys + use GFS_typedefs, only: GFS_control_type, GFS_grid_type, GFS_statein_type + use physcons, only: rhc_max + + type(GFS_control_type), intent(in) :: Model + type(GFS_grid_type), intent(in) :: Grid + type(GFS_statein_type), intent(in) :: Statein + + real(kind=kind_phys), intent(in) :: rhbbot, rhbtop + real(kind=kind_phys), dimension(size(Grid%xlon,1)), intent(in) :: work1, work2 + real(kind=kind_phys), intent(inout) :: clw(:,:,:), cnvc(:,:), cnvw(:,:) + integer, dimension(size(Grid%xlon,1)), intent(inout) :: ktop, kbot + real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levs), intent(out) :: rhc + + integer :: i,k + real(kind=kind_phys) :: tem + + clw(:,:,1) = 0.0 + clw(:,:,2) = -999.9 + if ((Model%imfdeepcnv >= 0) .or. (Model%imfshalcnv > 0)) then + cnvc(:,:) = 0.0 + cnvw(:,:) = 0.0 + endif + + ktop(:) = 1 + kbot(:) = Model%levs + + if (Model%ntcw > 0) then + do k=1,Model%levs + do i=1, size(Grid%xlon,1) + tem = rhbbot - (rhbbot-rhbtop) * (1.0-Statein%prslk(i,k)) + tem = rhc_max * work1(i) + tem * work2(i) + rhc(i,k) = max(0.0, min(1.0,tem)) + enddo + enddo + endif + +end subroutine GFS_suite_interstitial_4_run + +end module GFS_suite_interstitial_4 + +module GFS_suite_interstitial_5 + +contains + +subroutine GFS_suite_interstitial_5_init () +end subroutine GFS_suite_interstitial_5_init + +subroutine GFS_suite_interstitial_5_finalize() +end subroutine GFS_suite_interstitial_5_finalize + +!> \section arg_table_GFS_suite_interstitial_5_run Argument Table +!! | local var name | longname | description | units | rank | type | kind | intent | optional | +!! |----------------|--------------------------------------------------------|-----------------------------------------------------------------------|---------------|------|-------------------------------|-----------|--------|----------| +!! | clw | convective_transportable_tracers | array to contain cloud water and other convective trans. tracers | kg kg-1 | 3 | real | kind_phys | inout | F | +!! | cnvc | convective_cloud_cover | convective cloud cover | frac | 2 | real | kind_phys | inout | F | +!! | cnvw | convective_cloud_water_specific_humidity | convective cloud water specific humidity | kg kg-1 | 2 | real | kind_phys | inout | F | +!! +subroutine GFS_suite_interstitial_5_run (clw, cnvc, cnvw) + + use machine, only: kind_phys + + real(kind=kind_phys), allocatable, intent(inout) :: clw(:,:,:), cnvc(:,:), cnvw(:,:) + + deallocate (clw) + if (allocated(cnvc)) deallocate(cnvc) + if (allocated(cnvw)) deallocate(cnvw) + +end subroutine GFS_suite_interstitial_5_run + +end module diff --git a/physics/moninedmf.f b/physics/moninedmf.f index bf3314a25..55450e871 100755 --- a/physics/moninedmf.f +++ b/physics/moninedmf.f @@ -2,6 +2,10 @@ !! Contains most of the hybrid eddy-diffusivity mass-flux scheme except for the !! subroutine that calculates the mass flux and updraft properties. + module edmf + + contains + !> \defgroup HEDMF Hybrid Eddy-diffusivity Mass-flux PBL and Free Atmospheric Turbulence !! @{ !! \brief The Hybrid EDMF scheme is a first-order turbulent transport scheme used for subgrid-scale vertical turbulent mixing in the PBL and above. It blends the traditional first-order approach that has been used and improved over the last several years with a more recent scheme that uses a mass-flux approach to calculate the countergradient diffusion terms. @@ -13,67 +17,70 @@ !! \section intraphysics Intraphysics Communication !! This space is reserved for a description of how this scheme uses information from other scheme types and/or how information calculated in this scheme is used in other scheme types. + subroutine edmf_init () + end subroutine edmf_init + + subroutine edmf_finalize () + end subroutine edmf_finalize + !> \brief This subroutine contains all of logic for the Hybrid EDMF PBL scheme except for the calculation of the updraft properties and mass flux. !! !! The scheme works on a basic level by calculating background diffusion coefficients and updating them according to which processes are occurring in the column. The most important difference in diffusion coefficients occurs between those levels in the PBL and those above the PBL, so the PBL height calculation is of utmost importance. An initial estimate is calculated in a "predictor" step in order to calculate Monin-Obukhov similarity values and a corrector step recalculates the PBL height based on updated surface thermal characteristics. Using the PBL height and the similarity parameters, the diffusion coefficients are updated below the PBL top based on Hong and Pan (1996) \cite hong_and_pan_1996 (including counter-gradient terms). Diffusion coefficients in the free troposphere (above the PBL top) are calculated according to Louis (1979) \cite louis_1979 with updated Richardson number-dependent functions. If it is diagnosed that PBL top-down mixing is occurring according to Lock et al. (2000) \cite lock_et_al_2000 , then then diffusion coefficients are updated accordingly. Finally, for convective boundary layers (defined as when the Obukhov length exceeds a threshold), the counter-gradient terms are replaced using the mass flux scheme of Siebesma et al. (2007) \cite siebesma_et_al_2007 . In order to return time tendencies, a fully implicit solution is found using tridiagonal matrices, and time tendencies are "backed out." Before returning, the time tendency of temperature is updated to reflect heating due to TKE dissipation following Han et al. (2015) \cite han_et_al_2015 . !! -!! \section arg_table_edmf_run Arguments -!! | local var name | longname | description | units | rank | type | kind | intent | optional | -!! |----------------|-------------------------------------------------------|------------------------------------|---------|------|---------|-----------|--------|----------| -!! | im | horizontal_loop_extent | horizontal loop extent, start at 1 | index | 0 | integer | | in | F | -!! -!! \param[in] ix horizontal dimension -!! \param[in] im number of used points -!! \param[in] km vertical layer dimension -!! \param[in] ntrac number of tracers -!! \param[in] ntcw cloud condensate index in the tracer array -!! \param[in,out] dv v-momentum tendency (\f$ m s^{-2} \f$) -!! \param[in,out] du u-momentum tendency (\f$ m s^{-2} \f$) -!! \param[in,out] tau temperature tendency (\f$ K s^{-1} \f$) -!! \param[in,out] rtg moisture tendency (\f$ kg kg^{-1} s^{-1} \f$) -!! \param[in] u1 u component of layer wind (\f$ m s^{-1} \f$) -!! \param[in] v1 v component of layer wind (\f$ m s^{-1} \f$) -!! \param[in] t1 layer mean temperature (\f$ K \f$) -!! \param[in] q1 layer mean tracer concentration (units?) -!! \param[in] swh total sky shortwave heating rate (\f$ K s^-1 \f$) -!! \param[in] hlw total sky longwave heating rate (\f$ K s^-1 \f$) -!! \param[in] xmu time step zenith angle adjust factor for shortwave -!! \param[in] psk Exner function at surface interface? -!! \param[in] rbsoil surface bulk Richardson number -!! \param[in] zorl surface roughness (units?) -!! \param[in] u10m 10-m u wind (\f$ m s^{-1} \f$) -!! \param[in] v10m 10-m v wind (\f$ m s^{-1} \f$) -!! \param[in] fm fm parameter from PBL scheme -!! \param[in] fh fh parameter from PBL scheme -!! \param[in] tsea ground surface temperature (K) -!! \param[in] qss surface saturation humidity (units?) -!! \param[in] heat surface sensible heat flux (units?) -!! \param[in] evap evaporation from latent heat flux (units?) -!! \param[in] stress surface wind stress? (\f$ cm*v^2\f$ in sfc_diff subroutine) (units?) -!! \param[in] spd1 surface wind speed? (units?) -!! \param[out] kpbl PBL top index -!! \param[in] prsi pressure at layer interfaces (units?) -!! \param[in] del pressure difference between level k and k+1 (units?) -!! \param[in] prsl mean layer pressure (units?) -!! \param[in] prslk Exner function at layer -!! \param[in] phii interface geopotential height (units?) -!! \param[in] phil layer geopotential height (units?) -!! \param[in] delt physics time step (s) -!! \param[in] dspheat flag for TKE dissipative heating -!! \param[out] dusfc surface u-momentum tendency (units?) -!! \param[out] dvsfc surface v-momentum tendency (units?) -!! \param[out] dtsfc surface temperature tendency (units?) -!! \param[out] dqsfc surface moisture tendency (units?) -!! \param[out] hpbl PBL top height (m) -!! \param[out] hgamt counter gradient mixing term for temperature (units?) -!! \param[out] hgamq counter gradient mixing term for moisture (units?) -!! \param[out] dkt diffusion coefficient for temperature (units?) -!! \param[in] kinver index location of temperature inversion -!! \param[in] xkzm_m background vertical diffusion coefficient for momentum (units?) -!! \param[in] xkzm_h background vertical diffusion coefficeint for heat, moisture (units?) -!! \param[in] xkzm_s sigma threshold for background momentum diffusion (units?) -!! \param[in] lprnt flag to print some output -!! \param[in] ipr index of point to print +!! \section arg_table_edmf_run +!! | local var name | longname | description | units | rank | type | kind | intent | optional | +!! |----------------|--------------------------------------------------------------|----------------------------------------------------|---------------|------|---------|-----------|--------|----------| +!! | ix | horizontal_dimension | horizontal dimension | index | 0 | integer | | in | F | +!! | im | horizontal_loop_extent | horizontal loop extent, start at 1 | index | 0 | integer | | in | F | +!! | km | vertical_dimension | vertical layer dimension | index | 0 | integer | | in | F | +!! | ntrac | number_of_vertical_diffusion_tracers | number of tracers to diffuse vertically | count | 0 | integer | | in | F | +!! | ntcw | index_for_liquid_cloud_condensate | cloud condensate index in tracer array | index | 0 | integer | | in | F | +!! | dv | tendency_of_y_wind_due_to_model_physics | updated tendency of the y wind | m s-2 | 2 | real | kind_phys | inout | F | +!! | du | tendency_of_x_wind_due_to_model_physics | updated tendency of the x wind | m s-2 | 2 | real | kind_phys | inout | F | +!! | tau | tendency_of_air_temperature_due_to_model_physics | updated tendency of the temperature | K s-1 | 2 | real | kind_phys | inout | F | +!! | rtg | tendency_of_tracers_due_to_model_physics | updated tendency of the tracers | kg kg-1 s-1 | 3 | real | kind_phys | inout | F | +!! | u1 | x_wind | x component of layer wind | m s-1 | 2 | real | kind_phys | in | F | +!! | v1 | y_wind | y component of layer wind | m s-1 | 2 | real | kind_phys | in | F | +!! | t1 | air_temperature | layer mean air temperature | K | 2 | real | kind_phys | in | F | +!! | q1 | tracer_concentration | layer mean tracer concentration | kg kg-1 | 3 | real | kind_phys | in | F | +!! | swh | tendency_of_air_temperature_due_to_shortwave_heating | total sky shortwave heating rate | K s-1 | 2 | real | kind_phys | in | F | +!! | hlw | tendency_of_air_temperature_due_to_longwave_heating | total sky longwave heating rate | K s-1 | 2 | real | kind_phys | in | F | +!! | xmu | zenith_angle_temporal_adjustment_factor_for_shortwave_fluxes | zenith angle temporal adjustment factor for shortwave | none | 1 | real | kind_phys | in | F | +!! | psk | exner_function_at_lowest_model_interface | exner function at the surface interface | none | 1 | real | kind_phys | in | F | +!! | rbsoil | bulk_richardson_number_at_lowest_model_level | bulk Richardson number at the surface | none | 1 | real | kind_phys | in | F | +!! | zorl | surface_roughness_length | surface roughness length in cm | cm | 1 | real | kind_phys | in | F | +!! | u10m | x_wind_at_10m | x component of wind at 10 m | m s-1 | 1 | real | kind_phys | in | F | +!! | v10m | y_wind_at_10m | y component of wind at 10 m | m s-1 | 1 | real | kind_phys | in | F | +!! | fm | Monin-Obukhov_similarity_function_for_momentum | Monin-Obukhov similarity function for momentum | none | 1 | real | kind_phys | in | F | +!! | fh | Monin-Obukhov_similarity_function_for_heat | Monin-Obukhov similarity function for heat | none | 1 | real | kind_phys | in | F | +!! | tsea | surface_skin_temperature | surface temperature | K | 1 | real | kind_phys | in | F | +!! | heat | kinematic_surface_upward_sensible_heat_flux | kinematic surface upward sensible heat flux | K m s-1 | 1 | real | kind_phys | in | F | +!! | evap | kinematic_surface_upward_latent_heat_flux | kinematic surface upward latent heat flux | kg kg-1 m s-1 | 1 | real | kind_phys | in | F | +!! | stress | surface_wind_stress | surface wind stress | m2 s-2 | 1 | real | kind_phys | in | F | +!! | spd1 | wind_speed_at_lowest_model_layer | wind speed at lowest model level | m s-1 | 1 | real | kind_phys | in | F | +!! | kpbl | vertical_index_at_top_of_atmosphere_boundary_layer | PBL top model level index | index | 1 | integer | | out | F | +!! | prsi | air_pressure_at_interface | air pressure at model layer interfaces | Pa | 2 | real | kind_phys | in | F | +!! | del | air_pressure_difference_between_midlayers | pres(k) - pres(k+1) | Pa | 2 | real | kind_phys | in | F | +!! | prsl | air_pressure | mean layer pressure | Pa | 2 | real | kind_phys | in | F | +!! | prslk | dimensionless_exner_function | Exner function at layers | none | 2 | real | kind_phys | in | F | +!! | phii | geopotential_at_interface | geopotential at model layer interfaces | m2 s-2 | 2 | real | kind_phys | in | F | +!! | phil | geopotential | geopotential at model layer centers | m2 s-2 | 2 | real | kind_phys | in | F | +!! | delt | time_step_for_physics | time step for physics | s | 0 | real | kind_phys | in | F | +!! | dspheat | flag_TKE_dissipation_heating | flag for using TKE dissipation heating | flag | 0 | logical | | in | F | +!! | dusfc | instantaneous_surface_x_momentum_flux | x momentum flux | Pa | 1 | real | kind_phys | out | F | +!! | dvsfc | instantaneous_surface_y_momentum_flux | y momentum flux | Pa | 1 | real | kind_phys | out | F | +!! | dtsfc | instantaneous_surface_upward_sensible_heat_flux | surface upward sensible heat flux | W m-2 | 1 | real | kind_phys | out | F | +!! | dqsfc | instantaneous_surface_upward_latent_heat_flux | surface upward latent heat flux | W m-2 | 1 | real | kind_phys | out | F | +!! | hpbl | atmosphere_boundary_layer_thickness | PBL thickness | m | 1 | real | kind_phys | out | F | +!! | hgamt | countergradient_mixing_term_for_temperature | countergradient mixing term for temperature | K | 1 | real | kind_phys | out | F | +!! | hgamq | countergradient_mixing_term_for_water_vapor | countergradient mixing term for water vapor | kg kg-1 | 1 | real | kind_phys | out | F | +!! | dkt | atmosphere_heat_diffusivity | diffusivity for heat | m2 s-1 | 1 | real | kind_phys | out | F | +!! | kinver | index_of_highest_temperature_inversion | index of highest temperature inversion | index | 1 | integer | | in | F | +!! | xkzm_m | atmosphere_momentum_diffusivity_background | background value of momentum diffusivity | m2 s-1 | 0 | real | kind_phys | in | F | +!! | xkzm_h | atmosphere_heat_diffusivity_background | background value of heat diffusivity | m2 s-1 | 0 | real | kind_phys | in | F | +!! | xkzm_s | diffusivity_background_sigma_level | sigma level threshold for background diffusivity | none | 0 | real | kind_phys | in | F | +!! | lprnt | flag_print | flag for printing diagnostics to output | flag | 0 | logical | | in | F | +!! | ipr | horizontal_index_of_printed_column | horizontal index of printed column | index | 0 | integer | | in | F | !! !! \section general General Algorithm !! -# Compute preliminary variables from input arguments. @@ -91,10 +98,10 @@ !! -# Solve for the horizontal momentum tendencies and add them to output tendency terms. !! \section detailed Detailed Algorithm !! @{ - subroutine moninedmf(ix,im,km,ntrac,ntcw,dv,du,tau,rtg, & + subroutine edmf_run (ix,im,km,ntrac,ntcw,dv,du,tau,rtg, & & u1,v1,t1,q1,swh,hlw,xmu, & & psk,rbsoil,zorl,u10m,v10m,fm,fh, & - & tsea,qss,heat,evap,stress,spd1,kpbl, & + & tsea,heat,evap,stress,spd1,kpbl, & & prsi,del,prsl,prslk,phii,phil,delt,dspheat, & & dusfc,dvsfc,dtsfc,dqsfc,hpbl,hgamt,hgamq,dkt, & & kinver,xkzm_m,xkzm_h,xkzm_s,lprnt,ipr) @@ -121,7 +128,7 @@ subroutine moninedmf(ix,im,km,ntrac,ntcw,dv,du,tau,rtg, & & rbsoil(im), zorl(im), & & u10m(im), v10m(im), & & fm(im), fh(im), & - & tsea(im), qss(im), & + & tsea(im), & & spd1(im), & & prsi(ix,km+1), del(ix,km), & & prsl(ix,km), prslk(ix,km), & @@ -1039,7 +1046,7 @@ subroutine moninedmf(ix,im,km,ntrac,ntcw,dv,du,tau,rtg, & ! ! solve tridiagonal problem for heat and moisture ! -!> The tridiagonal system is solved by calling the internal ::tridin subroutine. +!> The tridiagonal system is solved by calling the internal ::edmf_tridin subroutine. call tridin(im,km,ntrac,al,ad,au,a1,a2,au,a1,a2) ! @@ -1189,7 +1196,7 @@ subroutine moninedmf(ix,im,km,ntrac,ntcw,dv,du,tau,rtg, & ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! return - end + end subroutine edmf_run !> @} c----------------------------------------------------------------------- @@ -1197,111 +1204,113 @@ subroutine moninedmf(ix,im,km,ntrac,ntcw,dv,du,tau,rtg, & !! \brief Routine to solve the tridiagonal system to calculate temperature and moisture at \f$ t + \Delta t \f$; part of two-part process to calculate time tendencies due to vertical diffusion. !! !! Origin of subroutine unknown. - subroutine tridi2(l,n,cl,cm,cu,r1,r2,au,a1,a2) -cc - use machine , only : kind_phys - implicit none - integer k,n,l,i - real(kind=kind_phys) fk -cc - real(kind=kind_phys) cl(l,2:n),cm(l,n),cu(l,n-1),r1(l,n),r2(l,n), & - & au(l,n-1),a1(l,n),a2(l,n) -c----------------------------------------------------------------------- - do i=1,l - fk = 1./cm(i,1) - au(i,1) = fk*cu(i,1) - a1(i,1) = fk*r1(i,1) - a2(i,1) = fk*r2(i,1) - enddo - do k=2,n-1 - do i=1,l - fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1)) - au(i,k) = fk*cu(i,k) - a1(i,k) = fk*(r1(i,k)-cl(i,k)*a1(i,k-1)) - a2(i,k) = fk*(r2(i,k)-cl(i,k)*a2(i,k-1)) - enddo - enddo - do i=1,l - fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1)) - a1(i,n) = fk*(r1(i,n)-cl(i,n)*a1(i,n-1)) - a2(i,n) = fk*(r2(i,n)-cl(i,n)*a2(i,n-1)) - enddo - do k=n-1,1,-1 - do i=1,l - a1(i,k) = a1(i,k)-au(i,k)*a1(i,k+1) - a2(i,k) = a2(i,k)-au(i,k)*a2(i,k+1) - enddo - enddo -c----------------------------------------------------------------------- - return - end -c----------------------------------------------------------------------- +C subroutine edmf_tridi2(l,n,cl,cm,cu,r1,r2,au,a1,a2) +c +C use machine , only : kind_phys +C implicit none +C integer k,n,l,i +C real(kind=kind_phys) fk +c +C real(kind=kind_phys) cl(l,2:n),cm(l,n),cu(l,n-1),r1(l,n),r2(l,n), & +C & au(l,n-1),a1(l,n),a2(l,n) +C----------------------------------------------------------------------- +C do i=1,l +C fk = 1./cm(i,1) +C au(i,1) = fk*cu(i,1) +C a1(i,1) = fk*r1(i,1) +C a2(i,1) = fk*r2(i,1) +C enddo +C do k=2,n-1 +C do i=1,l +C fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1)) +C au(i,k) = fk*cu(i,k) +C a1(i,k) = fk*(r1(i,k)-cl(i,k)*a1(i,k-1)) +C a2(i,k) = fk*(r2(i,k)-cl(i,k)*a2(i,k-1)) +C enddo +C enddo +C do i=1,l +C fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1)) +C a1(i,n) = fk*(r1(i,n)-cl(i,n)*a1(i,n-1)) +C a2(i,n) = fk*(r2(i,n)-cl(i,n)*a2(i,n-1)) +C enddo +C do k=n-1,1,-1 +C do i=1,l +C a1(i,k) = a1(i,k)-au(i,k)*a1(i,k+1) +C a2(i,k) = a2(i,k)-au(i,k)*a2(i,k+1) +C enddo +C enddo +C----------------------------------------------------------------------- +C return +C end +C +C----------------------------------------------------------------------- !> \ingroup HEDMF !! \brief Routine to solve the tridiagonal system to calculate u- and v-momentum at \f$ t + \Delta t \f$; part of two-part process to calculate time tendencies due to vertical diffusion. !! !! Origin of subroutine unknown. - subroutine tridin(l,n,nt,cl,cm,cu,r1,r2,au,a1,a2) -cc - use machine , only : kind_phys - implicit none - integer is,k,kk,n,nt,l,i - real(kind=kind_phys) fk(l) -cc - real(kind=kind_phys) cl(l,2:n), cm(l,n), cu(l,n-1), & - & r1(l,n), r2(l,n*nt), & - & au(l,n-1), a1(l,n), a2(l,n*nt), & - & fkk(l,2:n-1) -c----------------------------------------------------------------------- - do i=1,l - fk(i) = 1./cm(i,1) - au(i,1) = fk(i)*cu(i,1) - a1(i,1) = fk(i)*r1(i,1) - enddo - do k = 1, nt - is = (k-1) * n - do i = 1, l - a2(i,1+is) = fk(i) * r2(i,1+is) - enddo - enddo - do k=2,n-1 - do i=1,l - fkk(i,k) = 1./(cm(i,k)-cl(i,k)*au(i,k-1)) - au(i,k) = fkk(i,k)*cu(i,k) - a1(i,k) = fkk(i,k)*(r1(i,k)-cl(i,k)*a1(i,k-1)) - enddo - enddo - do kk = 1, nt - is = (kk-1) * n - do k=2,n-1 - do i=1,l - a2(i,k+is) = fkk(i,k)*(r2(i,k+is)-cl(i,k)*a2(i,k+is-1)) - enddo - enddo - enddo - do i=1,l - fk(i) = 1./(cm(i,n)-cl(i,n)*au(i,n-1)) - a1(i,n) = fk(i)*(r1(i,n)-cl(i,n)*a1(i,n-1)) - enddo - do k = 1, nt - is = (k-1) * n - do i = 1, l - a2(i,n+is) = fk(i)*(r2(i,n+is)-cl(i,n)*a2(i,n+is-1)) - enddo - enddo - do k=n-1,1,-1 - do i=1,l - a1(i,k) = a1(i,k) - au(i,k)*a1(i,k+1) - enddo - enddo - do kk = 1, nt - is = (kk-1) * n - do k=n-1,1,-1 - do i=1,l - a2(i,k+is) = a2(i,k+is) - au(i,k)*a2(i,k+is+1) - enddo - enddo - enddo -c----------------------------------------------------------------------- - return - end +C subroutine edmf_tridin(l,n,nt,cl,cm,cu,r1,r2,au,a1,a2) +c +C use machine , only : kind_phys +C implicit none +C integer is,k,kk,n,nt,l,i +C real(kind=kind_phys) fk(l) +c +C real(kind=kind_phys) cl(l,2:n), cm(l,n), cu(l,n-1), & +C & r1(l,n), r2(l,n*nt), & +C & au(l,n-1), a1(l,n), a2(l,n*nt), & +C & fkk(l,2:n-1) +C----------------------------------------------------------------------- +C do i=1,l +C fk(i) = 1./cm(i,1) +C au(i,1) = fk(i)*cu(i,1) +C a1(i,1) = fk(i)*r1(i,1) +C enddo +C do k = 1, nt +C is = (k-1) * n +C do i = 1, l +C a2(i,1+is) = fk(i) * r2(i,1+is) +C enddo +C enddo +C do k=2,n-1 +C do i=1,l +C fkk(i,k) = 1./(cm(i,k)-cl(i,k)*au(i,k-1)) +C au(i,k) = fkk(i,k)*cu(i,k) +C a1(i,k) = fkk(i,k)*(r1(i,k)-cl(i,k)*a1(i,k-1)) +C enddo +C enddo +C do kk = 1, nt +C is = (kk-1) * n +C do k=2,n-1 +C do i=1,l +C a2(i,k+is) = fkk(i,k)*(r2(i,k+is)-cl(i,k)*a2(i,k+is-1)) +C enddo +C enddo +C enddo +C do i=1,l +C fk(i) = 1./(cm(i,n)-cl(i,n)*au(i,n-1)) +C a1(i,n) = fk(i)*(r1(i,n)-cl(i,n)*a1(i,n-1)) +C enddo +C do k = 1, nt +C is = (k-1) * n +C do i = 1, l +C a2(i,n+is) = fk(i)*(r2(i,n+is)-cl(i,n)*a2(i,n+is-1)) +C enddo +C enddo +C do k=n-1,1,-1 +C do i=1,l +C a1(i,k) = a1(i,k) - au(i,k)*a1(i,k+1) +C enddo +C enddo +C do kk = 1, nt +C is = (kk-1) * n +C do k=n-1,1,-1 +C do i=1,l +C a2(i,k+is) = a2(i,k+is) - au(i,k)*a2(i,k+is+1) +C enddo +C enddo +C enddo +C----------------------------------------------------------------------- +C return +C end !> @} + end module edmf diff --git a/physics/tridi.f b/physics/tridi.f new file mode 100644 index 000000000..dbf03e62b --- /dev/null +++ b/physics/tridi.f @@ -0,0 +1,105 @@ + subroutine tridi2(l,n,cl,cm,cu,r1,r2,au,a1,a2) +cc + use machine , only : kind_phys + implicit none + integer k,n,l,i + real(kind=kind_phys) fk +cc + real(kind=kind_phys) cl(l,2:n),cm(l,n),cu(l,n-1),r1(l,n),r2(l,n), & + & au(l,n-1),a1(l,n),a2(l,n) +c---------------------------------------------------------------------- + do i=1,l + fk = 1./cm(i,1) + au(i,1) = fk*cu(i,1) + a1(i,1) = fk*r1(i,1) + a2(i,1) = fk*r2(i,1) + enddo + do k=2,n-1 + do i=1,l + fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1)) + au(i,k) = fk*cu(i,k) + a1(i,k) = fk*(r1(i,k)-cl(i,k)*a1(i,k-1)) + a2(i,k) = fk*(r2(i,k)-cl(i,k)*a2(i,k-1)) + enddo + enddo + do i=1,l + fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1)) + a1(i,n) = fk*(r1(i,n)-cl(i,n)*a1(i,n-1)) + a2(i,n) = fk*(r2(i,n)-cl(i,n)*a2(i,n-1)) + enddo + do k=n-1,1,-1 + do i=1,l + a1(i,k) = a1(i,k)-au(i,k)*a1(i,k+1) + a2(i,k) = a2(i,k)-au(i,k)*a2(i,k+1) + enddo + enddo +c----------------------------------------------------------------------- + return + end + +c----------------------------------------------------------------------- + + subroutine tridin(l,n,nt,cl,cm,cu,r1,r2,au,a1,a2) +cc + use machine , only : kind_phys + implicit none + integer is,k,kk,n,nt,l,i + real(kind=kind_phys) fk(l) +cc + real(kind=kind_phys) cl(l,2:n), cm(l,n), cu(l,n-1), & + & r1(l,n), r2(l,n*nt), & + & au(l,n-1), a1(l,n), a2(l,n*nt), & + & fkk(l,2:n-1) +c----------------------------------------------------------------------- + do i=1,l + fk(i) = 1./cm(i,1) + au(i,1) = fk(i)*cu(i,1) + a1(i,1) = fk(i)*r1(i,1) + enddo + do k = 1, nt + is = (k-1) * n + do i = 1, l + a2(i,1+is) = fk(i) * r2(i,1+is) + enddo + enddo + do k=2,n-1 + do i=1,l + fkk(i,k) = 1./(cm(i,k)-cl(i,k)*au(i,k-1)) + au(i,k) = fkk(i,k)*cu(i,k) + a1(i,k) = fkk(i,k)*(r1(i,k)-cl(i,k)*a1(i,k-1)) + enddo + enddo + do kk = 1, nt + is = (kk-1) * n + do k=2,n-1 + do i=1,l + a2(i,k+is) = fkk(i,k)*(r2(i,k+is)-cl(i,k)*a2(i,k+is-1)) + enddo + enddo + enddo + do i=1,l + fk(i) = 1./(cm(i,n)-cl(i,n)*au(i,n-1)) + a1(i,n) = fk(i)*(r1(i,n)-cl(i,n)*a1(i,n-1)) + enddo + do k = 1, nt + is = (k-1) * n + do i = 1, l + a2(i,n+is) = fk(i)*(r2(i,n+is)-cl(i,n)*a2(i,n+is-1)) + enddo + enddo + do k=n-1,1,-1 + do i=1,l + a1(i,k) = a1(i,k) - au(i,k)*a1(i,k+1) + enddo + enddo + do kk = 1, nt + is = (kk-1) * n + do k=n-1,1,-1 + do i=1,l + a2(i,k+is) = a2(i,k+is) - au(i,k)*a2(i,k+is+1) + enddo + enddo + enddo +c----------------------------------------------------------------------- + return + end