diff --git a/atmos_cubed_sphere/driver/fvGFS/atmosphere.F90 b/atmos_cubed_sphere/driver/fvGFS/atmosphere.F90 index a9cbaffbf..fe5aa5d66 100644 --- a/atmos_cubed_sphere/driver/fvGFS/atmosphere.F90 +++ b/atmos_cubed_sphere/driver/fvGFS/atmosphere.F90 @@ -1012,6 +1012,15 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block) enddo enddo enddo + if (.not. Atm(mytile)%flagstruct%hydrostatic) then + do k = 1, npz + do j = jsc,jec + do i = isc,iec + Atm(n)%delz(i,j,k) = Atm(n)%delz(i,j,k) + IAU_Data%delz_inc(i,j,k)*dt_atmos + enddo + enddo + enddo + endif ! add analysis increment to tracers to output from physics do nb = 1,Atm_block%nblks !if (nb.EQ.1) print*,'in block_update',IAU_Data%in_interval,IAU_Data%temp_inc(isc,jsc,30) diff --git a/atmos_cubed_sphere/tools/fv_iau_mod.F90 b/atmos_cubed_sphere/tools/fv_iau_mod.F90 index 23082b2a0..d63ef42d2 100644 --- a/atmos_cubed_sphere/tools/fv_iau_mod.F90 +++ b/atmos_cubed_sphere/tools/fv_iau_mod.F90 @@ -60,6 +60,7 @@ module fv_iau_mod real,allocatable :: va_inc(:,:,:) real,allocatable :: temp_inc(:,:,:) real,allocatable :: delp_inc(:,:,:) + real,allocatable :: delz_inc(:,:,:) real,allocatable :: tracer_inc(:,:,:,:) end type iau_internal_data_type type iau_external_data_type @@ -67,6 +68,7 @@ module fv_iau_mod real,allocatable :: va_inc(:,:,:) real,allocatable :: temp_inc(:,:,:) real,allocatable :: delp_inc(:,:,:) + real,allocatable :: delz_inc(:,:,:) real,allocatable :: tracer_inc(:,:,:,:) logical :: in_interval = .false. end type iau_external_data_type @@ -207,12 +209,14 @@ subroutine IAU_initialize (IPD_Control, IAU_Data,Init_parm) allocate(IAU_Data%va_inc(is:ie, js:je, km)) allocate(IAU_Data%temp_inc(is:ie, js:je, km)) allocate(IAU_Data%delp_inc(is:ie, js:je, km)) + allocate(IAU_Data%delz_inc(is:ie, js:je, km)) allocate(IAU_Data%tracer_inc(is:ie, js:je, km,ntracers)) ! allocate arrays that will hold iau state allocate (iau_state%inc1%ua_inc(is:ie, js:je, km)) allocate (iau_state%inc1%va_inc(is:ie, js:je, km)) allocate (iau_state%inc1%temp_inc (is:ie, js:je, km)) allocate (iau_state%inc1%delp_inc (is:ie, js:je, km)) + allocate (iau_state%inc1%delz_inc (is:ie, js:je, km)) allocate (iau_state%inc1%tracer_inc(is:ie, js:je, km,ntracers)) iau_state%hr1=IPD_Control%iaufhrs(1) call read_iau_forcing(IPD_Control,iau_state%inc1,'INPUT/'//trim(IPD_Control%iau_inc_files(1))) @@ -224,6 +228,7 @@ subroutine IAU_initialize (IPD_Control, IAU_Data,Init_parm) allocate (iau_state%inc2%va_inc(is:ie, js:je, km)) allocate (iau_state%inc2%temp_inc (is:ie, js:je, km)) allocate (iau_state%inc2%delp_inc (is:ie, js:je, km)) + allocate (iau_state%inc2%delz_inc (is:ie, js:je, km)) allocate (iau_state%inc2%tracer_inc(is:ie, js:je, km,ntracers)) iau_state%hr2=IPD_Control%iaufhrs(2) call read_iau_forcing(IPD_Control,iau_state%inc2,'INPUT/'//trim(IPD_Control%iau_inc_files(2))) @@ -302,6 +307,7 @@ subroutine updateiauforcing(IPD_Control,IAU_Data) IAU_Data%va_inc(i,j,k) =(delt*IAU_state%inc1%va_inc(i,j,k) + (1.-delt)* IAU_state%inc2%va_inc(i,j,k))*rdt IAU_Data%temp_inc(i,j,k) =(delt*IAU_state%inc1%temp_inc(i,j,k) + (1.-delt)* IAU_state%inc2%temp_inc(i,j,k))*rdt IAU_Data%delp_inc(i,j,k) =(delt*IAU_state%inc1%delp_inc(i,j,k) + (1.-delt)* IAU_state%inc2%delp_inc(i,j,k))*rdt + IAU_Data%delz_inc(i,j,k) =(delt*IAU_state%inc1%delz_inc(i,j,k) + (1.-delt)* IAU_state%inc2%delz_inc(i,j,k))*rdt do l=1,ntracers IAU_Data%tracer_inc(i,j,k,l) =(delt*IAU_state%inc1%tracer_inc(i,j,k,l) + (1.-delt)* IAU_state%inc2%tracer_inc(i,j,k,l))*rdt enddo @@ -327,6 +333,7 @@ subroutine setiauforcing(IPD_Control,IAU_Data) IAU_Data%va_inc(i,j,k) =IAU_state%inc1%va_inc(i,j,k)*rdt IAU_Data%temp_inc(i,j,k) =IAU_state%inc1%temp_inc(i,j,k)*rdt IAU_Data%delp_inc(i,j,k) =IAU_state%inc1%delp_inc(i,j,k)*rdt + IAU_Data%delz_inc(i,j,k) =IAU_state%inc1%delz_inc(i,j,k)*rdt do l = 1,ntracers IAU_Data%tracer_inc(i,j,k,l) =IAU_state%inc1%tracer_inc(i,j,k,l)*rdt enddo @@ -382,6 +389,7 @@ subroutine read_iau_forcing(IPD_Control,increments,fname) ! read in 1 time level call interp_inc('T_inc',increments%temp_inc(:,:,:),jbeg,jend) call interp_inc('delp_inc',increments%delp_inc(:,:,:),jbeg,jend) + call interp_inc('delz_inc',increments%delz_inc(:,:,:),jbeg,jend) call interp_inc('u_inc',increments%ua_inc(:,:,:),jbeg,jend) ! can these be treated as scalars? call interp_inc('v_inc',increments%va_inc(:,:,:),jbeg,jend) do l=1,ntracers diff --git a/atmos_cubed_sphere/tools/fv_treat_da_inc.F90 b/atmos_cubed_sphere/tools/fv_treat_da_inc.F90 index 0249423e5..a7200708b 100644 --- a/atmos_cubed_sphere/tools/fv_treat_da_inc.F90 +++ b/atmos_cubed_sphere/tools/fv_treat_da_inc.F90 @@ -172,6 +172,9 @@ subroutine read_da_inc(Atm, fv_domain) call apply_inc_on_3d_scalar('T_inc',Atm(1)%pt) call apply_inc_on_3d_scalar('delp_inc',Atm(1)%delp) + if (.not. Atm(1)%flagstruct%hydrostatic) then + call apply_inc_on_3d_scalar('delz_inc',Atm(1)%delz) + endif call apply_inc_on_3d_scalar('sphum_inc',Atm(1)%q(:,:,:,sphum)) call apply_inc_on_3d_scalar('liq_wat_inc',Atm(1)%q(:,:,:,liq_wat)) call apply_inc_on_3d_scalar('o3mr_inc',Atm(1)%q(:,:,:,o3mr)) diff --git a/gfsphysics/GFS_layer/GFS_diagnostics.F90 b/gfsphysics/GFS_layer/GFS_diagnostics.F90 index 161c592dd..498af7525 100644 --- a/gfsphysics/GFS_layer/GFS_diagnostics.F90 +++ b/gfsphysics/GFS_layer/GFS_diagnostics.F90 @@ -3354,6 +3354,24 @@ subroutine diag_populate (IPD_Diag, Model, Statein, Stateout, Sfcprop, Coupling, IPD_Diag(idx)%data(nb)%var3p => Diag(nb)%skebv_wts(:,:) enddo + !---dividing streamline + idx = idx + 1 + IPD_Diag(idx)%name = 'zmtnblck' + IPD_Diag(idx)%output_name = 'zmtnblck' + IPD_Diag(idx)%mod_name = 'physics' + IPD_Diag(idx)%file_name = 'surface' + IPD_Diag(idx)%desc = 'perturbation velocity' + IPD_Diag(idx)%unit = 'm/s' + IPD_Diag(idx)%type_stat_proc = 'inst' + IPD_Diag(idx)%level_type = '' + IPD_Diag(idx)%level = 1 + IPD_Diag(idx)%cnvfac = cn_one + IPD_Diag(idx)%zhour = Model%zhour + IPD_Diag(idx)%fcst_hour = Model%fhour + do nb = 1,nblks + IPD_Diag(idx)%data(nb)%var2p => Diag(nb)%zmtnblck(:) + enddo + !---sppt idx = idx + 1 IPD_Diag(idx)%name = 'sppt_wts' diff --git a/gfsphysics/GFS_layer/GFS_driver.F90 b/gfsphysics/GFS_layer/GFS_driver.F90 index 595bccc6e..a4404fa6e 100644 --- a/gfsphysics/GFS_layer/GFS_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_driver.F90 @@ -361,13 +361,13 @@ subroutine GFS_time_vary_step (Model, Statein, Stateout, Sfcprop, Coupling, & enddo enddo endif - if (Model%do_sppt) then - do nb = 1,nblks - do k=1,Model%levs - Diag(nb)%sppt_wts(:,k)=Coupling(nb)%sppt_wts(:,Model%levs-k+1) - enddo - enddo - endif + !if (Model%do_sppt) then + ! do nb = 1,nblks + ! do k=1,Model%levs + ! Diag(nb)%sppt_wts(:,k)=Coupling(nb)%sppt_wts(:,Model%levs-k+1) + ! enddo + ! enddo + !endif if (Model%do_shum) then do nb = 1,nblks do k=1,Model%levs @@ -407,14 +407,34 @@ subroutine GFS_stochastic_driver (Model, Statein, Stateout, Sfcprop, Coupling, & type(GFS_diag_type), intent(inout) :: Diag !--- local variables integer :: k, i - real(kind=kind_phys) :: upert, vpert, tpert, qpert, qnew + real(kind=kind_phys) :: upert, vpert, tpert, qpert, qnew,sppt_vwt if (Model%do_sppt) then do k = 1,size(Statein%tgrs,2) do i = 1,size(Statein%tgrs,1) + sppt_vwt=1.0 + if (Diag%zmtnblck(i).EQ.0.0) then + sppt_vwt=1.0 + else + if (k.GT.Diag%zmtnblck(i)+2) then + sppt_vwt=1.0 + endif + if (k.LE.Diag%zmtnblck(i)) then + sppt_vwt=0.0 + endif + if (k.EQ.Diag%zmtnblck(i)+1) then + sppt_vwt=0.333333 + endif + if (k.EQ.Diag%zmtnblck(i)+2) then + sppt_vwt=0.666667 + endif + endif + if (Model%use_zmtnblck)then + Coupling%sppt_wts(i,k)=(Coupling%sppt_wts(i,k)-1)*sppt_vwt+1.0 + endif + Diag%sppt_wts(i,Model%levs-k+1)=Coupling%sppt_wts(i,k) upert = (Stateout%gu0(i,k) - Statein%ugrs(i,k)) * Coupling%sppt_wts(i,k) vpert = (Stateout%gv0(i,k) - Statein%vgrs(i,k)) * Coupling%sppt_wts(i,k) -! tpert = (Stateout%gt0(i,k) - Statein%tgrs(i,k)) * Coupling%sppt_wts(i,k) - Tbd%dtdtr(i,k) tpert = (Stateout%gt0(i,k) - Statein%tgrs(i,k) - Tbd%dtdtr(i,k)) * Coupling%sppt_wts(i,k) qpert = (Stateout%gq0(i,k,1) - Statein%qgrs(i,k,1)) * Coupling%sppt_wts(i,k) diff --git a/gfsphysics/GFS_layer/GFS_physics_driver.F90 b/gfsphysics/GFS_layer/GFS_physics_driver.F90 index e3c1c3107..535970fd8 100644 --- a/gfsphysics/GFS_layer/GFS_physics_driver.F90 +++ b/gfsphysics/GFS_layer/GFS_physics_driver.F90 @@ -1510,7 +1510,8 @@ subroutine GFS_physics_driver & Sfcprop%hprime(1,1), oc, oa4, clx, theta, & sigma, gamma, elvmax, dusfcg, dvsfcg, & con_g, con_cp, con_rd, con_rv, Model%lonr, & - Model%nmtvr, Model%cdmbgwd, me, lprnt,ipr) + Model%nmtvr, Model%cdmbgwd, me, lprnt,ipr, & + Diag%zmtnblck) ! if (lprnt) print *,' dudtg=',dudt(ipr,:) diff --git a/gfsphysics/GFS_layer/GFS_typedefs.F90 b/gfsphysics/GFS_layer/GFS_typedefs.F90 index 9caa923ef..588883b1f 100644 --- a/gfsphysics/GFS_layer/GFS_typedefs.F90 +++ b/gfsphysics/GFS_layer/GFS_typedefs.F90 @@ -555,6 +555,7 @@ module GFS_typedefs !--- stochastic physics control parameters logical :: do_sppt + logical :: use_zmtnblck logical :: do_shum logical :: do_skeb integer :: skeb_npass @@ -842,10 +843,11 @@ module GFS_typedefs real (kind=kind_phys), pointer :: wet1 (:) => null() !< normalized soil wetness real (kind=kind_phys), pointer :: sr (:) => null() !< snow ratio : ratio of snow to total precipitation - real (kind=kind_phys), pointer :: skebu_wts(:,:) => null() !< 10 meater u/v wind speed - real (kind=kind_phys), pointer :: skebv_wts(:,:) => null() !< 10 meater u/v wind speed + real (kind=kind_phys), pointer :: skebu_wts(:,:) => null() !< 10 meater u/v wind speed + real (kind=kind_phys), pointer :: skebv_wts(:,:) => null() !< 10 meater u/v wind speed real (kind=kind_phys), pointer :: sppt_wts(:,:) => null() !< 10 meater u/v wind speed real (kind=kind_phys), pointer :: shum_wts(:,:) => null() !< 10 meater u/v wind speed + real (kind=kind_phys), pointer :: zmtnblck(:) => null() ! null() !< u momentum change due to physics real (kind=kind_phys), pointer :: dv3dt (:,:,:) => null() !< v momentum change due to physics @@ -1558,9 +1560,10 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & logical :: debug = .false. logical :: pre_rad = .false. !< flag for testing purpose !--- stochastic physics control parameters - logical :: do_sppt = .false. - logical :: do_shum = .false. - logical :: do_skeb = .false. + logical :: do_sppt = .false. + logical :: use_zmtnblck = .false. + logical :: do_shum = .false. + logical :: do_skeb = .false. integer :: skeb_npass = 11 !--- END NAMELIST VARIABLES @@ -1813,6 +1816,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & !--- stochastic physics options Model%do_sppt = do_sppt + Model%use_zmtnblck = use_zmtnblck Model%do_shum = do_shum Model%do_skeb = do_skeb @@ -2304,11 +2308,6 @@ subroutine control_print(Model) print *, ' xkzminv : ', Model%xkzminv print *, ' moninq_fac : ', Model%moninq_fac print *, ' ' - print *, 'stochastic physics' - print *, ' do_sppt : ', Model%do_sppt - print *, ' do_shum : ', Model%do_shum - print *, ' do_skeb : ', Model%do_skeb - print *, ' ' print *, 'tracers' print *, ' tracer_names : ', Model%tracer_names print *, ' ntrac : ', Model%ntrac @@ -2633,6 +2632,7 @@ subroutine diag_create (Diag, IM, Model) allocate (Diag%skebv_wts(IM,Model%levs)) allocate (Diag%sppt_wts(IM,Model%levs)) allocate (Diag%shum_wts(IM,Model%levs)) + allocate (Diag%zmtnblck(IM)) !--- 3D diagnostics if (Model%ldiag3d) then allocate (Diag%du3dt (IM,Model%levs,4)) diff --git a/gfsphysics/physics/gwdps.f b/gfsphysics/physics/gwdps.f index 00d3f6da1..2778c3697 100644 --- a/gfsphysics/physics/gwdps.f +++ b/gfsphysics/physics/gwdps.f @@ -184,7 +184,7 @@ SUBROUTINE GWDPS(IM,IX,IY,KM,A,B,C,U1,V1,T1,Q1,KPBL, & & PRSI,DEL,PRSL,PRSLK,PHII, PHIL,DELTIM,KDT, & & HPRIME,OC,OA4,CLX4,THETA,SIGMA,GAMMA,ELVMAX, & & DUSFC,DVSFC,G, CP, RD, RV, IMX, & - & nmtvr, cdmbgwd, me, lprnt, ipr) + & nmtvr, cdmbgwd, me, lprnt, ipr, rdxzb) ! ! ******************************************************************** ! -----> I M P L E M E N T A T I O N V E R S I O N <---------- @@ -286,7 +286,7 @@ SUBROUTINE GWDPS(IM,IX,IY,KM,A,B,C,U1,V1,T1,Q1,KPBL, & & U1(IX,KM), V1(IX,KM), T1(IX,KM), & & Q1(IX,KM), PRSI(IX,KM+1), DEL(IX,KM), & & PRSL(IX,KM), PRSLK(IX,KM), PHIL(IX,KM), & - & PHII(IX,KM+1) + & PHII(IX,KM+1),RDXZB(IY) real(kind=kind_phys) OC(IM), OA4(IY,4), CLX4(IY,4) & &, HPRIME(IM) ! for lm mtn blocking @@ -408,6 +408,7 @@ SUBROUTINE GWDPS(IM,IX,IY,KM,A,B,C,U1,V1,T1,Q1,KPBL, & ! IF ( NMTVR .eq. 14) then ! ---- for lm and gwd calculation points + RDXZB(:) = 0 ipt = 0 npt = 0 DO I = 1,IM @@ -585,7 +586,10 @@ SUBROUTINE GWDPS(IM,IX,IY,KM,A,B,C,U1,V1,T1,Q1,KPBL, & EK(I) = 0.5 * UP(I) * UP(I) ! --- Dividing Stream lime is found when PE =exceeds EK. - IF ( PE(I) .ge. EK(I) ) IDXZB(I) = K + IF ( PE(I) .ge. EK(I) ) THEN + IDXZB(I) = K + RDXZB(J) = real(K,kind=kind_phys) + ENDIF ! --- Then mtn blocked flow is between Zb=k(IDXZB(I)) and surface ! !> - The dividing streamline height (idxzb), of a subgrid scale @@ -709,6 +713,7 @@ SUBROUTINE GWDPS(IM,IX,IY,KM,A,B,C,U1,V1,T1,Q1,KPBL, & ! do i=1,npt IDXZB(i) = 0 + RDXZB(i) = 0. enddo ENDIF ! diff --git a/gfsphysics/physics/h2ointerp.f90 b/gfsphysics/physics/h2ointerp.f90 index 5548142a8..b9cbc11a3 100755 --- a/gfsphysics/physics/h2ointerp.f90 +++ b/gfsphysics/physics/h2ointerp.f90 @@ -126,7 +126,7 @@ subroutine h2ointerpol(me,npts,idate,fhour,jindx1,jindx2,h2oplout,ddy) integer idat(8),jdat(8) ! real(kind=kind_phys) ddy(npts) - real(kind=kind_phys) h2oplout(levh2o,npts,h2o_coeff) + real(kind=kind_phys) h2oplout(npts,levh2o,h2o_coeff) real(kind=kind_phys) rinc(5), rjday integer jdow, jdoy, jday real(4) rinc4(5) @@ -152,26 +152,23 @@ subroutine h2ointerpol(me,npts,idate,fhour,jindx1,jindx2,h2oplout,ddy) jday = 0 call w3doxdat(jdat,jdow,jdoy,jday) rjday = jdoy + jdat(5) / 24. - if (rjday < h2o_time(1)) rjday = rjday+365. + if (rjday < h2o_time(1)) rjday = rjday + 365. ! n2 = timeh2o + 1 - do j=1,timeh2o + do j=2,timeh2o if (rjday < h2o_time(j)) then n2 = j exit endif enddo n1 = n2 - 1 - if (n1 <= 0) n1 = n1 + timeh2o - if (n2 > timeh2o) n2 = n2 - timeh2o - ! ! if (me .eq. 0) print *,' n1=',n1,' n2=',n2,' rjday=',rjday ! &,'h2o_time=',h2o_time(n1),h2o_time(n2) ! - tx1 = (h2o_time(n2) - rjday) / (h2o_time(n2) - h2o_time(n1)) tx2 = 1.0 - tx1 + if (n2 > timeh2o) n2 = n2 - timeh2o ! do nc=1,h2o_coeff do l=1,levh2o @@ -179,8 +176,8 @@ subroutine h2ointerpol(me,npts,idate,fhour,jindx1,jindx2,h2oplout,ddy) j1 = jindx1(j) j2 = jindx2(j) tem = 1.0 - ddy(j) - h2oplout(j,l,nc) = & - tx1*(tem*h2oplin(j1,l,nc,n1)+ddy(j)*h2oplin(j2,l,nc,n1)) & + h2oplout(j,l,nc) = & + tx1*(tem*h2oplin(j1,l,nc,n1)+ddy(j)*h2oplin(j2,l,nc,n1)) & + tx2*(tem*h2oplin(j1,l,nc,n2)+ddy(j)*h2oplin(j2,l,nc,n2)) enddo enddo diff --git a/gfsphysics/physics/ozinterp.f90 b/gfsphysics/physics/ozinterp.f90 index f00927711..680d3d7ea 100644 --- a/gfsphysics/physics/ozinterp.f90 +++ b/gfsphysics/physics/ozinterp.f90 @@ -122,31 +122,30 @@ SUBROUTINE ozinterpol(me,npts,IDATE,FHOUR,jindx1,jindx2,ozplout,ddy) USE MACHINE, ONLY : kind_phys USE OZNE_DEF implicit none - integer iday,j,j1,j2,l,npts,nc,n1,n2 + integer iday,j,j1,j2,l,npts,nc,n1,n2 real(kind=kind_phys) fhour,tem, tx1, tx2 ! integer JINDX1(npts), JINDX2(npts) - integer me,idate(4) - integer IDAT(8),JDAT(8) + integer me, idate(4), IDAT(8),JDAT(8) ! real(kind=kind_phys) DDY(npts) real(kind=kind_phys) ozplout(npts,levozp,oz_coeff) real(kind=kind_phys) RINC(5), rjday integer jdow, jdoy, jday real(4) rinc4(5) - integer w3kindreal,w3kindint -! - IDAT=0 - IDAT(1)=IDATE(4) - IDAT(2)=IDATE(2) - IDAT(3)=IDATE(3) - IDAT(5)=IDATE(1) - RINC=0. - RINC(2)=FHOUR + integer w3kindreal, w3kindint +! + IDAT = 0 + IDAT(1) = IDATE(4) + IDAT(2) = IDATE(2) + IDAT(3) = IDATE(3) + IDAT(5) = IDATE(1) + RINC = 0. + RINC(2) = FHOUR call w3kind(w3kindreal,w3kindint) - if(w3kindreal==4) then - rinc4=rinc + if(w3kindreal == 4) then + rinc4 = rinc CALL W3MOVDAT(RINC4,IDAT,JDAT) else CALL W3MOVDAT(RINC,IDAT,JDAT) @@ -157,26 +156,25 @@ SUBROUTINE ozinterpol(me,npts,IDATE,FHOUR,jindx1,jindx2,ozplout,ddy) jday = 0 call w3doxdat(jdat,jdow,jdoy,jday) rjday = jdoy + jdat(5) / 24. - IF (RJDAY .LT. oz_time(1)) RJDAY = RJDAY+365. + IF (RJDAY < oz_time(1)) RJDAY = RJDAY + 365. ! n2 = timeoz + 1 - do j=1,timeoz - if (rjday .lt. oz_time(j)) then + do j=2,timeoz + if (rjday < oz_time(j)) then n2 = j exit endif enddo n1 = n2 - 1 - if (n1 <= 0) n1 = n1 + timeoz - if (n2 > timeoz) n2 = n2 - timeoz - ! -! if (me .eq. 0) print *,' n1=',n1,' n2=',n2,' rjday=',rjday +! if (me == 0) print *,' n1=',n1,' n2=',n2,' rjday=',rjday ! &,'oz_time=',oz_time(n1),oz_time(n2) ! tx1 = (oz_time(n2) - rjday) / (oz_time(n2) - oz_time(n1)) tx2 = 1.0 - tx1 + + if (n2 > timeoz) n2 = n2 - timeoz ! do nc=1,oz_coeff DO L=1,levozp @@ -184,9 +182,9 @@ SUBROUTINE ozinterpol(me,npts,IDATE,FHOUR,jindx1,jindx2,ozplout,ddy) J1 = JINDX1(J) J2 = JINDX2(J) TEM = 1.0 - DDY(J) - ozplout(j,L,nc) = & - tx1*(TEM*ozplin(J1,L,nc,n1)+DDY(J)*ozplin(J2,L,nc,n1)) & - + tx2*(TEM*ozplin(J1,L,nc,n2)+DDY(J)*ozplin(J2,L,nc,n2)) + ozplout(j,L,nc) = & + tx1*(TEM*ozplin(J1,L,nc,n1)+DDY(J)*ozplin(J2,L,nc,n1)) & + + tx2*(TEM*ozplin(J1,L,nc,n2)+DDY(J)*ozplin(J2,L,nc,n2)) ENDDO ENDDO enddo diff --git a/gfsphysics/physics/radlw_main.f b/gfsphysics/physics/radlw_main.f index 701505296..e60a54a82 100644 --- a/gfsphysics/physics/radlw_main.f +++ b/gfsphysics/physics/radlw_main.f @@ -4671,80 +4671,43 @@ subroutine taugb05 jplp = jpl + 1 jmo3p = jmo3 + 1 - if (specparm < 0.125 .and. specparm1 < 0.125) then + if (specparm < 0.125) then p0 = fs - f_one p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = fs1 - f_one - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 + 2 id210 = ind0 +11 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 + 2 - id211 = ind1 +11 - elseif (specparm > 0.875 .and. specparm1 > 0.875) then + elseif (specparm > 0.875) then p0 = -fs p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = -fs1 - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 + 1 id010 = ind0 +10 id100 = ind0 id110 = ind0 + 9 id200 = ind0 - 1 id210 = ind0 + 8 - - id001 = ind1 + 1 - id011 = ind1 +10 - id101 = ind1 - id111 = ind1 + 9 - id201 = ind1 - 1 - id211 = ind1 + 8 else fk00 = f_one - fs fk10 = fs fk20 = f_zero - fk01 = f_one - fs1 - fk11 = fs1 - fk21 = f_zero - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 id210 = ind0 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 - id211 = ind1 endif fac000 = fk00 * fac00(k) @@ -4754,6 +4717,45 @@ subroutine taugb05 fac110 = fk10 * fac10(k) fac210 = fk20 * fac10(k) + if (specparm1 < 0.125) then + p1 = fs1 - f_one + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + 2 + id211 = ind1 +11 + elseif (specparm1 > 0.875) then + p1 = -fs1 + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + 1 + id011 = ind1 +10 + id101 = ind1 + id111 = ind1 + 9 + id201 = ind1 - 1 + id211 = ind1 + 8 + else + fk01 = f_one - fs1 + fk11 = fs1 + fk21 = f_zero + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + id211 = ind1 + endif + fac001 = fk01 * fac01(k) fac101 = fk11 * fac01(k) fac201 = fk21 * fac01(k) @@ -5035,80 +5037,43 @@ subroutine taugb07 adjcolco2 = colamt(k,2) endif - if (specparm < 0.125 .and. specparm1 < 0.125) then + if (specparm < 0.125) then p0 = fs - f_one p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = fs1 - f_one - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 + 2 id210 = ind0 +11 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 + 2 - id211 = ind1 +11 - elseif (specparm > 0.875 .and. specparm1 > 0.875) then + elseif (specparm > 0.875) then p0 = -fs p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = -fs1 - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 + 1 id010 = ind0 +10 id100 = ind0 id110 = ind0 + 9 id200 = ind0 - 1 id210 = ind0 + 8 - - id001 = ind1 + 1 - id011 = ind1 +10 - id101 = ind1 - id111 = ind1 + 9 - id201 = ind1 - 1 - id211 = ind1 + 8 else fk00 = f_one - fs fk10 = fs fk20 = f_zero - fk01 = f_one - fs1 - fk11 = fs1 - fk21 = f_zero - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 id210 = ind0 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 - id211 = ind1 endif fac000 = fk00 * fac00(k) @@ -5118,6 +5083,45 @@ subroutine taugb07 fac110 = fk10 * fac10(k) fac210 = fk20 * fac10(k) + if (specparm1 < 0.125) then + p1 = fs1 - f_one + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + 2 + id211 = ind1 +11 + elseif (specparm1 > 0.875) then + p1 = -fs1 + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + 1 + id011 = ind1 +10 + id101 = ind1 + id111 = ind1 + 9 + id201 = ind1 - 1 + id211 = ind1 + 8 + else + fk01 = f_one - fs1 + fk11 = fs1 + fk21 = f_zero + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + id211 = ind1 + endif + fac001 = fk01 * fac01(k) fac101 = fk11 * fac01(k) fac201 = fk21 * fac01(k) @@ -5419,81 +5423,43 @@ subroutine taugb09 adjcoln2o = colamt(k,4) endif - if (specparm < 0.125 .and. specparm1 < 0.125) then + if (specparm < 0.125) then p0 = fs - f_one p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = fs1 - f_one - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 + 2 id210 = ind0 +11 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 + 2 - id211 = ind1 +11 - - elseif (specparm > 0.875 .and. specparm1 > 0.875) then + elseif (specparm > 0.875) then p0 = -fs p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = -fs1 - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 + 1 id010 = ind0 +10 id100 = ind0 id110 = ind0 + 9 id200 = ind0 - 1 id210 = ind0 + 8 - - id001 = ind1 + 1 - id011 = ind1 +10 - id101 = ind1 - id111 = ind1 + 9 - id201 = ind1 - 1 - id211 = ind1 + 8 else fk00 = f_one - fs fk10 = fs fk20 = f_zero - fk01 = f_one - fs1 - fk11 = fs1 - fk21 = f_zero - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 id210 = ind0 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 - id211 = ind1 endif fac000 = fk00 * fac00(k) @@ -5503,6 +5469,45 @@ subroutine taugb09 fac110 = fk10 * fac10(k) fac210 = fk20 * fac10(k) + if (specparm1 < 0.125) then + p1 = fs1 - f_one + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + 2 + id211 = ind1 +11 + elseif (specparm1 > 0.875) then + p1 = -fs1 + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + 1 + id011 = ind1 +10 + id101 = ind1 + id111 = ind1 + 9 + id201 = ind1 - 1 + id211 = ind1 + 8 + else + fk01 = f_one - fs1 + fk11 = fs1 + fk21 = f_zero + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + id211 = ind1 + endif + fac001 = fk01 * fac01(k) fac101 = fk11 * fac01(k) fac201 = fk21 * fac01(k) @@ -5806,80 +5811,43 @@ subroutine taugb12 indfp = indf + 1 jplp = jpl + 1 - if (specparm < 0.125 .and. specparm1 < 0.125) then + if (specparm < 0.125) then p0 = fs - f_one p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = fs1 - f_one - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 + 2 id210 = ind0 +11 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 + 2 - id211 = ind1 +11 - elseif (specparm > 0.875 .and. specparm1 > 0.875) then + elseif (specparm > 0.875) then p0 = -fs p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = -fs1 - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 + 1 id010 = ind0 +10 id100 = ind0 id110 = ind0 + 9 id200 = ind0 - 1 id210 = ind0 + 8 - - id001 = ind1 + 1 - id011 = ind1 +10 - id101 = ind1 - id111 = ind1 + 9 - id201 = ind1 - 1 - id211 = ind1 + 8 else fk00 = f_one - fs fk10 = fs fk20 = f_zero - fk01 = f_one - fs1 - fk11 = fs1 - fk21 = f_zero - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 id210 = ind0 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 - id211 = ind1 endif fac000 = fk00 * fac00(k) @@ -5889,6 +5857,45 @@ subroutine taugb12 fac110 = fk10 * fac10(k) fac210 = fk20 * fac10(k) + if (specparm1 < 0.125) then + p1 = fs1 - f_one + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + 2 + id211 = ind1 +11 + elseif (specparm1 > 0.875) then + p1 = -fs1 + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + 1 + id011 = ind1 +10 + id101 = ind1 + id111 = ind1 + 9 + id201 = ind1 - 1 + id211 = ind1 + 8 + else + fk01 = f_one - fs1 + fk11 = fs1 + fk21 = f_zero + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + id211 = ind1 + endif + fac001 = fk01 * fac01(k) fac101 = fk11 * fac01(k) fac201 = fk21 * fac01(k) @@ -6035,80 +6042,43 @@ subroutine taugb13 adjcolco2 = colamt(k,2) endif - if (specparm < 0.125 .and. specparm1 < 0.125) then + if (specparm < 0.125) then p0 = fs - f_one p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = fs1 - f_one - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 + 2 id210 = ind0 +11 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 + 2 - id211 = ind1 +11 - elseif (specparm > 0.875 .and. specparm1 > 0.875) then + elseif (specparm > 0.875) then p0 = -fs p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = -fs1 - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 + 1 id010 = ind0 +10 id100 = ind0 id110 = ind0 + 9 id200 = ind0 - 1 id210 = ind0 + 8 - - id001 = ind1 + 1 - id011 = ind1 +10 - id101 = ind1 - id111 = ind1 + 9 - id201 = ind1 - 1 - id211 = ind1 + 8 else fk00 = f_one - fs fk10 = fs fk20 = f_zero - fk01 = f_one - fs1 - fk11 = fs1 - fk21 = f_zero - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 id210 = ind0 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 - id211 = ind1 endif fac000 = fk00 * fac00(k) @@ -6118,6 +6088,45 @@ subroutine taugb13 fac110 = fk10 * fac10(k) fac210 = fk20 * fac10(k) + if (specparm1 < 0.125) then + p1 = fs1 - f_one + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + 2 + id211 = ind1 +11 + elseif (specparm1 > 0.875) then + p1 = -fs1 + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + 1 + id011 = ind1 +10 + id101 = ind1 + id111 = ind1 + 9 + id201 = ind1 - 1 + id211 = ind1 + 8 + else + fk01 = f_one - fs1 + fk11 = fs1 + fk21 = f_zero + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + id211 = ind1 + endif + fac001 = fk01 * fac01(k) fac101 = fk11 * fac01(k) fac201 = fk21 * fac01(k) @@ -6325,82 +6334,44 @@ subroutine taugb15 indmp = indm + 1 jplp = jpl + 1 jmn2p = jmn2 + 1 - - if (specparm < 0.125 .and. specparm1 < 0.125) then + if (specparm < 0.125) then p0 = fs - f_one p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = fs1 - f_one - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 + 2 id210 = ind0 +11 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 + 2 - id211 = ind1 +11 - elseif (specparm > 0.875 .and. specparm1 > 0.875) then + elseif (specparm > 0.875) then p0 = -fs p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = -fs1 - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 + 1 id010 = ind0 +10 id100 = ind0 id110 = ind0 + 9 id200 = ind0 - 1 id210 = ind0 + 8 - - id001 = ind1 + 1 - id011 = ind1 +10 - id101 = ind1 - id111 = ind1 + 9 - id201 = ind1 - 1 - id211 = ind1 + 8 else fk00 = f_one - fs fk10 = fs fk20 = f_zero - fk01 = f_one - fs1 - fk11 = fs1 - fk21 = f_zero - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 id210 = ind0 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 - id211 = ind1 endif fac000 = fk00 * fac00(k) @@ -6410,6 +6381,45 @@ subroutine taugb15 fac110 = fk10 * fac10(k) fac210 = fk20 * fac10(k) + if (specparm1 < 0.125) then + p1 = fs1 - f_one + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + 2 + id211 = ind1 +11 + elseif (specparm1 > 0.875) then + p1 = -fs1 + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + 1 + id011 = ind1 +10 + id101 = ind1 + id111 = ind1 + 9 + id201 = ind1 - 1 + id211 = ind1 + 8 + else + fk01 = f_one - fs1 + fk11 = fs1 + fk21 = f_zero + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + id211 = ind1 + endif + fac001 = fk01 * fac01(k) fac101 = fk11 * fac01(k) fac201 = fk21 * fac01(k) @@ -6517,80 +6527,43 @@ subroutine taugb16 indfp = indf + 1 jplp = jpl + 1 - if (specparm < 0.125 .and. specparm1 < 0.125) then + if (specparm < 0.125) then p0 = fs - f_one p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = fs1 - f_one - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 + 2 id210 = ind0 +11 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 + 2 - id211 = ind1 +11 - elseif (specparm > 0.875 .and. specparm1 > 0.875) then + elseif (specparm > 0.875) then p0 = -fs p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = -fs1 - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 + 1 id010 = ind0 +10 id100 = ind0 id110 = ind0 + 9 id200 = ind0 - 1 id210 = ind0 + 8 - - id001 = ind1 + 1 - id011 = ind1 +10 - id101 = ind1 - id111 = ind1 + 9 - id201 = ind1 - 1 - id211 = ind1 + 8 else fk00 = f_one - fs fk10 = fs fk20 = f_zero - fk01 = f_one - fs1 - fk11 = fs1 - fk21 = f_zero - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 id210 = ind0 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 - id211 = ind1 endif fac000 = fk00 * fac00(k) @@ -6600,6 +6573,45 @@ subroutine taugb16 fac110 = fk10 * fac10(k) fac210 = fk20 * fac10(k) + if (specparm1 < 0.125) then + p1 = fs1 - f_one + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + 2 + id211 = ind1 +11 + elseif (specparm1 > 0.875) then + p1 = -fs1 + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + 1 + id011 = ind1 +10 + id101 = ind1 + id111 = ind1 + 9 + id201 = ind1 - 1 + id211 = ind1 + 8 + else + fk01 = f_one - fs1 + fk11 = fs1 + fk21 = f_zero + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + id211 = ind1 + endif + fac001 = fk01 * fac01(k) fac101 = fk11 * fac01(k) fac201 = fk21 * fac01(k) diff --git a/io/FV3GFS_io.F90 b/io/FV3GFS_io.F90 index b6c33deb7..b3a3231de 100644 --- a/io/FV3GFS_io.F90 +++ b/io/FV3GFS_io.F90 @@ -108,7 +108,7 @@ module FV3GFS_io_mod real(kind=kind_phys),dimension(:,:),allocatable :: lon real(kind=kind_phys),dimension(:,:),allocatable :: lat real(kind=kind_phys),dimension(:,:),allocatable :: uwork - logical :: uwork_set + logical :: uwork_set = .false. character(128) :: uwindname integer, parameter :: DIAG_SIZE = 500 ! real(kind=kind_phys), parameter :: missing_value = 1.d30 @@ -2167,6 +2167,7 @@ subroutine gfdl_diag_register(Time, Sfcprop, Cldprop, Gfs_diag, Gfs_grid, Atm_bl Diag(idx)%desc = 'averaged potential evaporation rate' Diag(idx)%unit = 'W/M**2' Diag(idx)%mod_name = 'gfs_phys' + Diag(idx)%time_avg = .TRUE. allocate (Diag(idx)%data(nblks)) do nb = 1,nblks Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%ep(:) @@ -2747,6 +2748,17 @@ subroutine gfdl_diag_register(Time, Sfcprop, Cldprop, Gfs_diag, Gfs_grid, Atm_bl Diag(idx)%data(nb)%var3 => Gfs_diag(nb)%skebv_wts(:,:) enddo + idx = idx + 1 + Diag(idx)%axes = 2 + Diag(idx)%name = 'zmtnblck' + Diag(idx)%desc = 'level of dividing streamline' + Diag(idx)%unit = 'm/s' + Diag(idx)%mod_name = 'gfs_phys' + allocate (Diag(idx)%data(nblks)) + do nb = 1,nblks + Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%zmtnblck(:) + enddo + idx = idx + 1 Diag(idx)%axes = 3 Diag(idx)%name = 'sppt_wts' @@ -2770,7 +2782,6 @@ subroutine gfdl_diag_register(Time, Sfcprop, Cldprop, Gfs_diag, Gfs_grid, Atm_bl enddo ! if(mpp_pe()==mpp_root_pe())print *,'in gfdl_diag_register,af shum_wts,idx=',idx - !--- three-dimensional variables that need to be handled special when writing !rab do num = 1,6 !rab write (xtra,'(I1)') num diff --git a/module_fcst_grid_comp.F90 b/module_fcst_grid_comp.F90 index 74ffab1ef..0bb48c81f 100644 --- a/module_fcst_grid_comp.F90 +++ b/module_fcst_grid_comp.F90 @@ -390,7 +390,7 @@ subroutine fcst_initialize(fcst_comp, importState, exportState, clock, rc) write(dateSH,'(I2.2)')date_init(4) write(dateSN,'(I2.2)')date_init(5) write(dateSS,'(I2.2)')date_init(6) - dateS="hours since "//dateSY//'-'//dateSM//'-'//dateSD//'-'//dateSH//':'// & + dateS="hours since "//dateSY//'-'//dateSM//'-'//dateSD//' '//dateSH//':'// & dateSN//":"//dateSS if(mype==0) print *,'dateS=',trim(dateS),'date_init=',date_init call ESMF_AttributeSet(exportState, convention="NetCDF", purpose="FV3", & diff --git a/stochastic_physics/compns_stochy.F90 b/stochastic_physics/compns_stochy.F90 index 1d6f03d3f..1cdeccd4a 100644 --- a/stochastic_physics/compns_stochy.F90 +++ b/stochastic_physics/compns_stochy.F90 @@ -48,7 +48,7 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) shum_lscale,fhstoch,stochini,skeb_varspect_opt,sppt_sfclimit, & skeb,skeb_tau,skeb_vdof,skeb_lscale,iseed_skeb,skeb_vfilt,skeb_diss_smooth, & skeb_sigtop1,skeb_sigtop2,skebnorm,sppt_sigtop1,sppt_sigtop2,& - shum_sigefold,skebint,skeb_npass + shum_sigefold,skebint,skeb_npass,use_zmtnblck tol=0.01 ! tolerance for calculations ! spectral resolution defintion ntrunc=-999 @@ -61,6 +61,7 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) skeb = -999. ! stochastic KE backscatter amplitude ! logicals do_sppt = .false. + use_zmtnblck = .false. do_shum = .false. do_skeb = .false. ! for SKEB random patterns. @@ -124,13 +125,16 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) IF (skeb(1) > 0 ) THEN do_skeb=.true. if (skebnorm==0) then ! stream function norm - skeb=skeb*5.0e5/sqrt(deltim) + skeb=skeb*1.111e3*sqrt(deltim) + !skeb=skeb*5.0e5/sqrt(deltim) endif if (skebnorm==1) then ! stream function norm - skeb=skeb*1/sqrt(deltim) + skeb=skeb*0.00222*sqrt(deltim) + !skeb=skeb*1/sqrt(deltim) endif if (skebnorm==2) then ! vorticty function norm - skeb=skeb*5.0e-7/sqrt(deltim) + skeb=skeb*1.111e-9*sqrt(deltim) + !skeb=skeb*5.0e-7/sqrt(deltim) endif ! adjust skeb values for resolution. ! scaling is such that a value of 1.0 at T574 with a 900 second diff --git a/stochastic_physics/stochastic_physics.F90 b/stochastic_physics/stochastic_physics.F90 index 10b81a1eb..b2d6f3d4c 100644 --- a/stochastic_physics/stochastic_physics.F90 +++ b/stochastic_physics/stochastic_physics.F90 @@ -32,6 +32,7 @@ subroutine init_stochastic_physics(Model,Init_parm) call init_stochdata(Model%levs,Model%dtp,Model%input_nml_file,Model%fn_nml,Init_parm%nlunit) ! check to see decomposition Model%do_sppt=do_sppt +Model%use_zmtnblck=use_zmtnblck Model%do_shum=do_shum Model%do_skeb=do_skeb Model%skeb_npass=skeb_npass @@ -54,7 +55,7 @@ subroutine init_stochastic_physics(Model,Init_parm) enddo if (sppt_sfclimit) then vfact_sppt(2)=vfact_sppt(3)*0.5 - vfact_sppt(1)=vfact_sppt(3)*0.25 + vfact_sppt(1)=0.0 endif if (is_master()) then do k=1,MOdel%levs diff --git a/stochastic_physics/stochy_namelist_def.F90 b/stochastic_physics/stochy_namelist_def.F90 index 3b84fd302..5a09fd58b 100644 --- a/stochastic_physics/stochy_namelist_def.F90 +++ b/stochastic_physics/stochy_namelist_def.F90 @@ -23,6 +23,6 @@ module stochy_namelist_def integer,dimension(5) ::skeb_vfilt integer(8),dimension(5) ::iseed_sppt,iseed_shum,iseed_skeb logical stochini,sppt_logit - logical do_shum,do_sppt,do_skeb + logical do_shum,do_sppt,do_skeb,use_zmtnblck end module stochy_namelist_def diff --git a/stochastic_physics/stochy_patterngenerator.F90 b/stochastic_physics/stochy_patterngenerator.F90 index 7c40e53c0..d155c7df2 100644 --- a/stochastic_physics/stochy_patterngenerator.F90 +++ b/stochastic_physics/stochy_patterngenerator.F90 @@ -42,7 +42,7 @@ subroutine patterngenerator_init(lscale, delt, tscale, stdev, iseed, rpattern,& integer, intent(in) :: nlon,nlat,jcap,npatterns,varspect_opt integer, intent(in) :: ls_node(ls_dim,3),nlevs type(random_pattern), intent(out), dimension(npatterns) :: rpattern - integer(8), intent(in) :: iseed(npatterns) + integer(8), intent(inout) :: iseed(npatterns) integer m,j,l,n,nm,nn,np,indev1,indev2,indod1,indod2 integer(8) count, count_rate, count_max, count_trunc integer(8) :: iscale = 10000000000 @@ -55,6 +55,15 @@ subroutine patterngenerator_init(lscale, delt, tscale, stdev, iseed, rpattern,& nlats = nlat ntrunc = jcap ndimspec = (ntrunc+1)*(ntrunc+2)/2 +! propagate seed supplied from namelist to all patterns... + if (iseed(1) .NE. 0) then + do np=2,npatterns + if (iseed(np).EQ.0) then + iseed(np)=iseed(1)+np*100000000 + endif + enddo + endif + do np=1,npatterns allocate(rpattern(np)%idx(0:ntrunc,0:ntrunc)) allocate(rpattern(np)%idx_e(len_trie_ls))