module tfilt_massfix contains subroutine tfilt_massfixrun (ztodt, lat, u3m1, u3, & v3m1, v3, t3m1, t3, q3m1, & q3, psm1, ps, alpha, & etamid, qfcst, vort, div, vortm2, & divm2, qminus, psm2, um2, & vm2, tm2, qm2, vortm1, divm1, & omga, dpsl, dpsm, beta, hadv , & nlon, pdeldry, pdelm1dry, pdelm2dry) !----------------------------------------------------------------------- ! ! Purpose: ! Time filter (second half of filter for vorticity and divergence only) ! ! Method: ! ! Author: ! !----------------------------------------------------------------------- use shr_kind_mod, only: r8 => shr_kind_r8 use cam_history, only: outfld use eul_control_mod, only : fixmas,eps use pmgrid, only: plon, plev, plevp, plat use pspect use commap use constituents, only: pcnst, qmin, cnst_cam_outfld, & tottnam, tendnam, cnst_get_type_byind, fixcnam, & hadvnam, vadvnam use time_manager, only: get_nstep use physconst, only: cpair, gravit use scamMod use phys_control, only: phys_getopts #if ( defined BFB_CAM_SCAM_IOP ) use iop use constituents, only: cnst_get_ind, cnst_name #endif implicit none ! ! Input arguments ! real(r8), intent(in) :: ztodt ! two delta t (unless nstep<2) real(r8), intent(inout) :: qfcst(plon,plev,pcnst)! slt moisture forecast real(r8), intent(in) :: vort(plon,plev) real(r8), intent(in) :: div(plon,plev) real(r8), intent(inout) :: vortm2(plon,plev) real(r8), intent(inout) :: divm2(plon,plev) real(r8), intent(in) :: qminus(plon,plev,pcnst) real(r8), intent(inout) :: psm2(plon) real(r8), intent(inout) :: um2(plon,plev) real(r8), intent(inout) :: vm2(plon,plev) real(r8), intent(inout) :: tm2(plon,plev) real(r8), intent(inout) :: qm2(plon,plev,pcnst) real(r8), intent(inout) :: omga(plon,plev) real(r8), intent(in) :: dpsl(plon) real(r8), intent(in) :: dpsm(plon) real(r8), intent(in) :: beta ! energy fixer coefficient real(r8), intent(in) :: hadv(plon,plev,pcnst) ! horizonal q advection tendency real(r8), intent(in) :: alpha(pcnst) real(r8), intent(in) :: etamid(plev) ! vertical coords at midpoints real(r8), intent(in) :: u3(plon,plev) real(r8), intent(in) :: v3(plon,plev) real(r8), intent(inout) :: t3(plon,plev) real(r8), intent(inout) :: pdeldry(:,:) ! dry pressure difference at time n3 real(r8), intent(inout) :: pdelm1dry(:,:) ! dry pressure difference at time n3m1 real(r8), intent(in) :: pdelm2dry(:,:) ! dry pressure difference at time n3m2 integer, intent(in) :: lat integer, intent(in) :: nlon ! Input/Output arguments real(r8), intent(inout) :: q3(plon,plev,pcnst) real(r8), intent(inout) :: ps(plon) real(r8), intent(inout) :: vortm1(plon,plev) real(r8), intent(inout) :: psm1(plon) real(r8), intent(inout) :: u3m1(plon,plev) real(r8), intent(inout) :: v3m1(plon,plev) real(r8), intent(inout) :: t3m1(plon,plev) real(r8), intent(inout) :: divm1(plon,plev) real(r8), intent(inout) :: q3m1(plon,plev,pcnst) ! ! Local workspace ! integer ifcnt ! Counter integer :: nstep ! current timestep number integer :: timefiltstep ! real(r8) tfix (plon) ! T correction real(r8) engycorr(plon,plev) ! energy equivalent to T correction real(r8) rpmid(plon,plev) ! 1./pmid real(r8) pdel(plon,plev) ! pdel(k) = pint (k+1)-pint (k) real(r8) pint(plon,plevp) ! pressure at model interfaces (n ) real(r8) pmid(plon,plev) ! pressure at model levels (time n) real(r8) utend(plon,plev) ! du/dt real(r8) vtend(plon,plev) ! dv/dt real(r8) ttend(plon,plev) ! dT/dt real(r8) qtend(plon,plev,pcnst)! dq/dt real(r8) pstend(plon) ! d(ps)/dt real(r8) vadv (plon,plev,pcnst) ! vertical q advection tendency real(r8) pintm1(plon,plevp) ! pressure at model interfaces (n-1) real(r8) pmidm1(plon,plev) ! pressure at model levels (time n-1) real(r8) pdelm1(plon,plev) ! pdelm1(k) = pintm1(k+1)-pintm1(k) real(r8) om2eps real(r8) corm real(r8) wm real(r8) absf real(r8) worst logical lfixlim ! flag to turn on fixer limiter real(r8) ta(plon,plev,pcnst) ! total advection of constituents real(r8) dqfx3(plon,plev,pcnst)! q tendency due to mass adjustment real(r8) coslat ! cosine(latitude) real(r8) rcoslat(plon) ! 1./cosine(latitude) ! real(r8) engt ! Thermal energy integral ! real(r8) engk ! Kinetic energy integral ! real(r8) engp ! Potential energy integral integer i, k, m,j,ixcldliq,ixcldice,ixnumliq,ixnumice #if ( defined BFB_CAM_SCAM_IOP ) real(r8) :: t3forecast(plon,plev),delta_t3(plon,plev) real(r8) :: q3forecast(plon,plev,pcnst),delta_q3(plon,plev,pcnst) #endif real(r8) fixmas_plon(plon) real(r8) beta_plon(plon) real(r8) clat_plon(plon) real(r8) alpha_plon(plon) character(len=16) :: microp_scheme ! microphysics scheme !----------------------------------------------------------------------- nstep = get_nstep() #if ( defined BFB_CAM_SCAM_IOP ) ! ! Calculate 3d dynamics term ! do k=1,plev do i=1,nlon divt3dsav(i,k,lat)=(t3(i,k)-tm2(i,k))/ztodt -t2sav(i,k,lat) t3forecast(i,k)=tm2(i,k)+ztodt*t2sav(i,k,lat)+ztodt*divt3dsav(i,k,lat) end do end do do i=1,nlon do m=1,pcnst do k=1,plev divq3dsav(i,k,m,lat)= (qfcst(i,k,m)-qminus(i,k,m))/ztodt q3forecast(i,k,m)=qminus(i,k,m)+divq3dsav(i,k,m,lat)*ztodt end do end do end do q3(:nlon,:,:)=q3forecast(:nlon,:,:) t3(:nlon,:)=t3forecast(:nlon,:) qfcst(:nlon,:,:)=q3(:nlon,:,:) ! ! outflds for iop history tape - to get bit for bit with scam ! the n-1 values are put out. After the fields are written out ! the current time level of info will be buffered for output next ! timestep ! call outfld('t',t3 ,plon ,lat ) call outfld('q',q3 ,plon ,lat ) call outfld('Ps',ps ,plon ,lat ) call outfld('u',u3 ,plon ,lat ) call outfld('v',v3 ,plon ,lat ) ! ! read single values into plon arrays for output to history tape ! it would be nice if history tape supported 1 dimensional array variables ! fixmas_plon(:)=fixmas beta_plon(:)=beta clat_plon(:)=clat(lat) call outfld('fixmas',fixmas_plon,plon ,lat ) call outfld('beta',beta_plon ,plon ,lat ) call outfld('CLAT ',clat_plon ,plon ,lat ) call outfld('divT3d',divt3dsav(1,1,lat) ,plon ,lat ) do m =1,pcnst call outfld(trim(cnst_name(m))//'_dten',divq3dsav(1,1,m,lat) ,plon ,lat ) end do #endif coslat = cos(clat(lat)) do i=1,nlon rcoslat(i) = 1._r8/coslat enddo lfixlim = .true. ! ! Set average dry mass to specified constant preserving horizontal ! gradients of ln(ps). Proportionality factor was calculated in STEPON ! for nstep=0 or SCAN2 otherwise from integrals calculated in INIDAT ! and SCAN2 respectively. ! Set p*. ! do i=1,nlon ps(i) = ps(i)*fixmas end do ! ! Set current time pressure arrays for model levels etc. ! call plevs0(nlon ,plon ,plev ,ps ,pint ,pmid ,pdel) ! !$OMP PARALLEL DO PRIVATE (K, I) do k=1,plev do i=1,nlon rpmid(i,k) = 1._r8/pmid(i,k) enddo enddo ! ! Add temperature correction for energy conservation ! !$OMP PARALLEL DO PRIVATE (K, I) do k=1,plev do i=1,nlon engycorr(i,k) = (cpair/gravit)*beta*pdel(i,k)/ztodt t3 (i,k) = t3(i,k) + beta end do end do do i=1,nlon tfix(i) = beta/ztodt end do ! ! Output Energy correction term ! ! using do loop and select in order to enable functional parallelism with OpenMP !$OMP PARALLEL DO PRIVATE (I) do i=1,2 select case (i) case (1) call outfld ('ENGYCORR',engycorr ,plon ,lat ) case (2) call outfld ('TFIX ',tfix ,plon ,lat ) end select end do ! ! Compute q tendency due to mass adjustment ! If LFIXLIM = .T., then: ! Check to see if fixer is exceeding a desired fractional limit of the ! constituent mixing ratio ("corm"). If so, then limit the fixer to ! that specified limit. ! do m=1,pcnst if (cnst_get_type_byind(m).eq.'dry' ) then corm = 1.e36_r8 else corm = 0.1_r8 end if !$OMP PARALLEL DO PRIVATE (K, I, IFCNT, WORST, WM, ABSF) do k=1,plev do i=1,nlon if (single_column) then dqfx3(i,k,m) = dqfxcam(i,k,m) else dqfx3(i,k,m) = alpha(m)*etamid(k)*abs(qfcst(i,k,m) - qminus(i,k,m)) #if ( defined BFB_CAM_SCAM_IOP ) dqfx3sav(i,k,m,lat) = dqfx3(i,k,m) #endif endif end do if (lfixlim) then ifcnt = 0 worst = 0._r8 wm = 0._r8 do i = 1,nlon absf = abs(dqfx3(i,k,m)) if (absf.gt.corm) then ifcnt = ifcnt + 1 worst = max(absf,worst) wm = wm + absf dqfx3(i,k,m) = sign(corm,dqfx3(i,k,m)) endif end do if (ifcnt.gt.0) then wm = wm/real(ifcnt,r8) ! TBH: Commented out as of CAM CRB meeting on 6/20/03 ! write(iulog,1000) m,corm,ifcnt,k,lat,wm,worst endif endif do i=1,nlon dqfx3(i,k,m) = qfcst(i,k,m)*dqfx3(i,k,m)/ztodt q3 (i,k,m) = qfcst(i,k,m) + ztodt*dqfx3(i,k,m) ta (i,k,m) = (q3 (i,k,m) - qminus(i,k,m))/ztodt vadv (i,k,m) = (qfcst(i,k,m) - qminus(i,k,m))/ztodt - hadv(i,k,m) end do end do end do !$OMP PARALLEL DO PRIVATE (K, I) do k=1,plev do i=1,nlon pdeldry(i,k) = pdel(i,k)*(1._r8-q3(i,k,1)) end do ! i end do ! k #if ( defined BFB_CAM_SCAM_IOP ) do m=1,pcnst alpha_plon(:)= alpha(m) call outfld(trim(cnst_name(m))//'_alph',alpha_plon ,plon ,lat ) call outfld(trim(cnst_name(m))//'_dqfx',dqfx3sav(1,1,m,lat) ,plon ,lat ) end do #endif ! ! Check for and correct invalid constituents ! call qneg3 ('TFILT_MASSFIX',lat ,nlon ,plon ,plev , & 1, pcnst, qmin ,q3(1,1,1)) ! ! Send slt tendencies to the history tape ! !$OMP PARALLEL DO PRIVATE (M) do m=1,pcnst if ( cnst_cam_outfld(m) ) then call outfld(tottnam(m),ta(1,1,m),plon ,lat ) end if end do if (.not. single_column) then ! ! Calculate vertical motion field ! call omcalc (rcoslat ,div ,u3 ,v3 ,dpsl , & dpsm ,pmid ,pdel ,rpmid ,pint(1,plevp), & omga ,nlon ) endif ! write(iulog,*)'tfilt: lat=',lat ! write(iulog,*)'omga=',omga ! ! Time filter (second half of filter for vorticity and divergence only) ! ! if(lat.eq.2) then ! write(iulog,*)'tfilt: ps=',psm2(13),psm1(13),ps(13) ! write(iulog,*)'tfilt: u=',um2(13,18),u3m1(13,18),u3(13,18) ! write(iulog,*)'tfilt: t=',tm2(13,18),t3m1(13,18),t3(13,18) ! write(iulog,*)'tfilt: water=',qm2(13,18,1),q3m1(13,18,1),q3(13,18,1) ! write(iulog,*)'tfilt: cwat=',qm2(13,18,2),q3m1(13,18,2),q3(13,18,2) ! write(iulog,*)'tfilt: vort=',vortm2(13,18),vortm1(13,18),vort(13,18) ! write(iulog,*)'tfilt: div=',divm2(13,18),divm1(13,18),div(13,18) ! end if om2eps = 1._r8 - 2._r8*eps if (nstep.ge.2) then !$OMP PARALLEL DO PRIVATE (K, I, M) do k=1,plev do i=1,nlon u3m1(i,k) = om2eps*u3m1(i,k) + eps*um2(i,k) + eps*u3(i,k) v3m1(i,k) = om2eps*v3m1(i,k) + eps*vm2(i,k) + eps*v3(i,k) t3m1(i,k) = om2eps*t3m1(i,k) + eps*tm2(i,k) + eps*t3(i,k) q3m1(i,k,1) = om2eps*q3m1(i,k,1) + eps*qm2(i,k,1) + eps*q3(i,k,1) vortm1(i,k) = om2eps*vortm1(i,k) + eps*vortm2(i,k) + eps*vort(i,k) divm1(i,k) = om2eps*divm1(i,k) + eps*divm2(i,k) + eps*div(i,k) end do do m=2,pcnst if (cnst_get_type_byind(m) .eq. 'wet') then do i=1,nlon q3m1(i,k,m) = om2eps*q3m1(i,k,m) + eps*qm2(i,k,m) + eps*q3(i,k,m) end do endif end do do m=2,pcnst if (cnst_get_type_byind(m) .eq. 'dry') then do i=1,nlon ! calculate numerator (timefiltered mass * pdeldry) q3m1(i,k,m) = (om2eps*pdelm1dry(i,k)*q3m1(i,k,m) + & eps*pdelm2dry(i,k)*qm2(i,k,m) + & eps*pdeldry(i,k)*q3(i,k,m)) end do !i endif !dry end do !m do i=1,nlon ! calculate time filtered value of pdeldry pdelm1dry(i,k) = om2eps*pdelm1dry(i,k) + & eps*pdelm2dry(i,k) + eps*pdeldry(i,k) end do !i ! divide time filtered mass*pdeldry by timefiltered pdeldry do m=2,pcnst if (cnst_get_type_byind(m) == 'dry') then do i=1,nlon q3m1(i,k,m) = q3m1(i,k,m)/pdelm1dry(i,k) end do !i endif ! dry end do !m end do do i=1,nlon psm1(i) = om2eps*psm1(i) + eps*psm2(i) + eps*ps(i) end do end if call plevs0 (nlon ,plon ,plev ,psm1 ,pintm1 ,pmidm1 ,pdelm1) ! ! Compute time tendencies:comment out since currently not on h-t ! !$OMP PARALLEL DO PRIVATE (K, I) do k=1,plev do i=1,nlon ttend(i,k) = (t3(i,k)-tm2(i,k))/ztodt utend(i,k) = (u3(i,k)-um2(i,k))/ztodt vtend(i,k) = (v3(i,k)-vm2(i,k))/ztodt end do end do !$OMP PARALLEL DO PRIVATE (M, K, I) do m=1,pcnst do k=1,plev do i=1,nlon qtend(i,k,m) = (q3(i,k,m) - qm2(i,k,m))/ztodt end do end do end do do i=1,nlon pstend(i) = (ps(i) - psm2(i))/ztodt end do !$OMP PARALLEL DO PRIVATE (M) do m=1,pcnst if ( cnst_cam_outfld(m) ) then call outfld (tendnam(m),qtend(1,1,m),plon,lat) call outfld (fixcnam(m),dqfx3(1,1,m),plon,lat) call outfld (hadvnam(m),hadv (1,1,m),plon,lat) call outfld (vadvnam(m),vadv (1,1,m),plon,lat) end if end do ! using do loop and select in order to enable functional parallelism with OpenMP !$OMP PARALLEL DO PRIVATE (I) do i=1,4 select case (i) case (1) call outfld ('UTEND ',utend,plon,lat) case (2) call outfld ('VTEND ',vtend,plon,lat) case (3) call outfld ('TTEND ',ttend,plon,lat) case (4) call outfld ('LPSTEN ',pstend,plon,lat) end select end do return 1000 format(' TIMEFILTER: WARNING: fixer for tracer ',i3,' exceeded ', & f8.5,' for ',i5,' points at k,lat = ',2i4, & ' Avg/Worst = ',1p2e10.2) end subroutine tfilt_massfixrun end module tfilt_massfix