module convadj_mod implicit none contains subroutine convadj(ngrid,nlay,nq,ptimestep, & pplay,pplev,ppopsk, & pu,pv,ph,pq, & pdufi,pdvfi,pdhfi,pdqfi, & pduadj,pdvadj,pdhadj,pdqadj) use tracer_h, only: igcm_h2o_vap use comcstfi_mod, only: g use generic_cloud_common_h, only: epsi_generic use generic_tracer_index_mod, only: generic_tracer_index use callkeys_mod, only: tracer,water,generic_condensation, & virtual_correction implicit none !================================================================== ! ! Purpose ! ------- ! Compute dry convective adjustment. ! See old reference paper: Hourdin et al. JAS 1993 ! "Meteorological Variability and the Annual Surface ! Pressure Cycle on Mars" ! https://doi.org/10.1175/1520-0469(1993)050%3C3625:MVATAS%3E2.0.CO;2 !================================================================== ! Arguments ! --------- INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns INTEGER,INTENT(IN) :: nlay ! number of atmospheric layers INTEGER,INTENT(IN) :: nq ! number of tracers REAL,INTENT(IN) :: ptimestep ! physics time step (s) REAL,INTENT(IN) :: pplay(ngrid,nlay) ! mi-layer pressure (Pa) REAL,INTENT(IN) :: pplev(ngrid,nlay+1) ! inter-layer pressure (Pa) REAL,INTENT(IN) :: ppopsk(ngrid,nlay) ! Exner (wrt surface pressure) ! fields from dynamics REAL,INTENT(IN) :: ph(ngrid,nlay) ! potential temperature (K) REAL,INTENT(IN) :: pu(ngrid,nlay) ! zonal wind (m/s) REAL,INTENT(IN) :: pv(ngrid,nlay) ! meridioanl wind (m/s) REAL,INTENT(IN) :: pq(ngrid,nlay,nq) ! tracers (../kg_air) ! tendencies due to previous physical processes REAL,INTENT(IN) :: pdufi(ngrid,nlay) ! on zonal wind (m/s/s) REAL,INTENT(IN) :: pdvfi(ngrid,nlay) ! on meridional wind (m/s/s) REAL,INTENT(IN) :: pdhfi(ngrid,nlay)! on potential temperature (/K/s) REAL,INTENT(IN) :: pdqfi(ngrid,nlay,nq) ! on tracers (../kg_air/s) ! tendencies due to convetive adjustement REAL,INTENT(OUT) :: pduadj(ngrid,nlay) ! on zonal wind (m/s/s) REAL,INTENT(OUT) :: pdvadj(ngrid,nlay) ! on meridinal wind (m/s/s) REAL,INTENT(OUT) :: pdhadj(ngrid,nlay) ! on potential temperature (/K/s) REAL,INTENT(OUT) :: pdqadj(ngrid,nlay,nq) ! on traceurs (../kg_air/s) ! Local ! ----- INTEGER ig,i,l,l1,l2,jj INTEGER jcnt, jadrs(ngrid) REAL sig(nlay+1),sdsig(nlay),dsig(nlay) REAL zu(ngrid,nlay),zv(ngrid,nlay) REAL zh(ngrid,nlay), zvh(ngrid,nlay) REAL zu2(ngrid,nlay),zv2(ngrid,nlay) REAL zh2(ngrid,nlay),zvh2(ngrid,nlay),zhc(ngrid,nlay) REAL zhm,zsm,zdsm,zum,zvm,zalpha,zhmc ! Tracers INTEGER iq REAL zq(ngrid,nlay,nq), zq2(ngrid,nlay,nq) REAL zqm(nq) integer igcm_generic_vap, igcm_generic_ice logical call_ice_vap_generic LOGICAL vtest(ngrid),down ! for conservation test real masse,cadjncons ! -------------- ! Initialisation ! -------------- zh(:,:)=ph(:,:)+pdhfi(:,:)*ptimestep zu(:,:)=pu(:,:)+pdufi(:,:)*ptimestep zv(:,:)=pv(:,:)+pdvfi(:,:)*ptimestep if(tracer) then zq(:,:,:)=pq(:,:,:)+pdqfi(:,:,:)*ptimestep end if zh2(:,:)=zh(:,:) zu2(:,:)=zu(:,:) zv2(:,:)=zv(:,:) zq2(:,:,:)=zq(:,:,:) ! ----------------------------- ! Detection of unstable columns ! ----------------------------- ! If ph(above) < ph(below) we set vtest=.true. DO ig=1,ngrid vtest(ig)=.false. ENDDO if((generic_condensation) .and. (virtual_correction)) THEN DO iq=1,nq call generic_tracer_index(nq,iq,igcm_generic_vap, & igcm_generic_ice,call_ice_vap_generic) if(call_ice_vap_generic) then zvh(:,:)=zh(:,:)* & (1.+zq(:,:,igcm_generic_vap)/epsi_generic)/ & (1.+zq(:,:,igcm_generic_vap)) endif ENDDO zvh2(:,:)=zvh(:,:) zhc(:,:)=zvh2(:,:) else zhc(:,:)=zh2(:,:) endif ! Find out which grid points are convectively unstable DO l=2,nlay DO ig=1,ngrid IF (zhc(ig,l).LT.zhc(ig,l-1)) THEN vtest(ig)=.true. ENDIF ENDDO ENDDO ! Make a list of them jcnt=0 DO ig=1,ngrid IF(vtest(ig)) THEN jcnt=jcnt+1 jadrs(jcnt)=ig ENDIF ENDDO ! --------------------------------------------------------------- ! Adjustment of the "jcnt" unstable profiles indicated by "jadrs" ! --------------------------------------------------------------- DO jj = 1, jcnt ! loop on every convective grid point i = jadrs(jj) ! Calculate sigma in this column DO l=1,nlay+1 sig(l)=pplev(i,l)/pplev(i,1) ENDDO DO l=1,nlay dsig(l)=sig(l)-sig(l+1) sdsig(l)=ppopsk(i,l)*dsig(l) ENDDO l2 = 1 ! Test loop upwards on l2 DO l2 = l2 + 1 IF (l2 .GT. nlay) EXIT IF (zhc(i, l2).LT.zhc(i, l2-1)) THEN ! l2 is the highest level of the unstable column l1 = l2 - 1 l = l1 zsm = sdsig(l2) zdsm = dsig(l2) zhm = zh2(i, l2) ! Test loop downwards DO zsm = zsm + sdsig(l) zdsm = zdsm + dsig(l) zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm if (generic_condensation .and. virtual_correction) then zhmc = zhm* & (1.+zq(i,l,igcm_generic_vap)/epsi_generic)/ & (1.+zq(i,l,igcm_generic_vap)) else zhmc = zhm endif ! do we have to extend the column downwards? down = .false. IF (l1 .ne. 1) then !-- and then IF (zhmc.LT.zhc(i, l1-1)) then down = .true. END IF END IF ! this could be a problem... if (down) then l1 = l1 - 1 l = l1 else ! can we extend the column upwards? if (l2 .eq. nlay) exit if (zhc(i, l2+1) .ge. zhmc) exit l2 = l2 + 1 l = l2 end if enddo ! New constant profile (average value) zalpha=0. zum=0. zvm=0. do iq=1,nq zqm(iq) = 0. end do DO l = l1, l2 zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l) zh2(i, l) = zhm ! modifs by RDW !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! zum=zum+dsig(l)*zu2(i,l) zvm=zvm+dsig(l)*zv2(i,l) ! zum=zum+dsig(l)*zu(i,l) ! zvm=zvm+dsig(l)*zv(i,l) do iq=1,nq zqm(iq) = zqm(iq)+dsig(l)*zq2(i,l,iq) ! zqm(iq) = zqm(iq)+dsig(l)*zq(i,l,iq) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! to conserve tracers/ KE, we must calculate zum, zvm and zqm using ! the up-to-date column values. If we do not do this, there are cases ! where convection stops at one level and starts at the next where we ! can break conservation of stuff (particularly tracers) significantly. end do ENDDO zalpha=zalpha/(zhm*(sig(l1)-sig(l2+1))) zum=zum/(sig(l1)-sig(l2+1)) zvm=zvm/(sig(l1)-sig(l2+1)) do iq=1,nq zqm(iq) = zqm(iq)/(sig(l1)-sig(l2+1)) end do IF(zalpha.GT.1.) THEN zalpha=1. ELSE ! IF(zalpha.LT.0.) STOP IF(zalpha.LT.1.e-4) zalpha=1.e-4 ENDIF DO l=l1,l2 zu2(i,l)=zu2(i,l)+zalpha*(zum-zu2(i,l)) zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l)) do iq=1,nq ! zq2(i,l,iq)=zq2(i,l,iq)+zalpha*(zqm(iq)-zq2(i,l,iq)) zq2(i,l,iq)=zqm(iq) end do ENDDO l2 = l2 + 1 END IF ! End of l1 to l2 instability treatment ! We now continue to test from l2 upwards ENDDO ! End of upwards loop on l2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! check conservation cadjncons=0.0 if(water)then do l = 1, nlay masse = (pplev(i,l) - pplev(i,l+1))/g iq = igcm_h2o_vap cadjncons = cadjncons + & masse*(zq2(i,l,iq)-zq(i,l,iq))/ptimestep end do endif if(cadjncons.lt.-1.e-6)then print*,'convadj has just crashed...' print*,'i = ',i print*,'l1 = ',l1 print*,'l2 = ',l2 print*,'cadjncons = ',cadjncons do l = 1, nlay print*,'dsig = ',dsig(l) end do do l = 1, nlay print*,'dsig = ',dsig(l) end do do l = 1, nlay print*,'sig = ',sig(l) end do do l = 1, nlay print*,'pplay(ig,:) = ',pplay(i,l) end do do l = 1, nlay+1 print*,'pplev(ig,:) = ',pplev(i,l) end do do l = 1, nlay print*,'ph(ig,:) = ',ph(i,l) end do do l = 1, nlay print*,'ph(ig,:) = ',ph(i,l) end do do l = 1, nlay print*,'ph(ig,:) = ',ph(i,l) end do do l = 1, nlay print*,'zh(ig,:) = ',zh(i,l) end do do l = 1, nlay print*,'zh2(ig,:) = ',zh2(i,l) end do do l = 1, nlay print*,'zq(ig,:,vap) = ',zq(i,l,igcm_h2o_vap) end do do l = 1, nlay print*,'zq2(ig,:,vap) = ',zq2(i,l,igcm_h2o_vap) end do print*,'zqm(vap) = ',zqm(igcm_h2o_vap) print*,'jadrs=',jadrs call abort endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ENDDO pdhadj(:,:)=(zh2(:,:)-zh(:,:))/ptimestep pduadj(:,:)=(zu2(:,:)-zu(:,:))/ptimestep pdvadj(:,:)=(zv2(:,:)-zv(:,:))/ptimestep if(tracer) then pdqadj(:,:,:)=(zq2(:,:,:)-zq(:,:,:))/ptimestep end if end subroutine convadj end module convadj_mod