subroutine condens_n2(ngrid,nlayer,nq,ptimestep, & pcapcal,pplay,pplev,ptsrf,pt, & pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq, & picen2,psolaralb,pemisurf, & pdtc,pdtsrfc,pdpsrf,pduc,pdvc, & pdqc,pdicen2) use radinc_h, only : naerkind use comgeomfi_h implicit none !================================================================== ! Purpose ! ------- ! Condense and/or sublime N2 ice on the ground and in the ! atmosphere, and sediment the ice. ! ! Inputs ! ------ ! ngrid Number of vertical columns ! nlayer Number of layers ! pplay(ngrid,nlayer) Pressure layers ! pplev(ngrid,nlayer+1) Pressure levels ! pt(ngrid,nlayer) Temperature (in K) ! ptsrf(ngrid) Surface temperature ! ! pdt(ngrid,nlayermx) Time derivative before condensation/sublimation of pt ! pdtsrf(ngrid) Time derivative before condensation/sublimation of ptsrf ! picen2(ngrid) n2 ice at the surface (kg/m2) ! ! Outputs ! ------- ! pdpsrf(ngrid) \ Contribution of condensation/sublimation ! pdtc(ngrid,nlayermx) / to the time derivatives of Ps, pt, and ptsrf ! pdtsrfc(ngrid) / ! pdicen2(ngrid) Tendancy n2 ice at the surface (kg/m2) ! ! Both ! ---- ! psolaralb(ngrid) Albedo at the surface ! pemisurf(ngrid) Emissivity of the surface ! ! Authors ! ------- ! Francois Forget (1996,2013) ! Converted to Fortran 90 and slightly modified by R. Wordsworth (2009) ! ! !================================================================== #include "dimensions.h" #include "dimphys.h" #include "comcstfi.h" #include "surfdat.h" #include "comvert.h" #include "callkeys.h" #include "tracer.h" !----------------------------------------------------------------------- ! Arguments INTEGER ngrid, nlayer, nq REAL ptimestep REAL pplay(ngrid,nlayer),pplev(ngrid,nlayer+1) REAL pcapcal(ngrid) REAL pt(ngrid,nlayer) REAL ptsrf(ngrid),flu1(ngrid),flu2(ngrid),flu3(ngrid) REAL pphi(ngrid,nlayer) REAL pdt(ngrid,nlayer),pdtsrf(ngrid),pdtc(ngrid,nlayer) REAL pdtsrfc(ngrid),pdpsrf(ngrid) REAL picen2(ngrid),psolaralb(ngrid),pemisurf(ngrid) REAL pu(ngrid,nlayer) , pv(ngrid,nlayer) REAL pdu(ngrid,nlayer) , pdv(ngrid,nlayer) REAL pduc(ngrid,nlayer) , pdvc(ngrid,nlayer) REAL pq(ngridmx,nlayer,nq),pdq(ngrid,nlayer,nq) REAL pdqc(ngrid,nlayer,nq) !----------------------------------------------------------------------- ! Local variables INTEGER l,ig,ilay,it,iq,ich4_gas REAL*8 zt(ngridmx,nlayermx) REAL tcond_n2 REAL pcond_n2 REAL glob_average2d ! function REAL zqn2(ngridmx,nlayermx) ! N2 MMR used to compute Tcond/zqn2 REAL ztcond (ngridmx,nlayermx) REAL ztcondsol(ngridmx),zfallheat REAL pdicen2(ngridmx) REAL zcondicea(ngridmx,nlayermx), zcondices(ngridmx) REAL zfallice(ngridmx,nlayermx+1) REAL zmflux(nlayermx+1) REAL zu(nlayermx),zv(nlayermx) REAL zq(nlayermx,nqmx),zq1(nlayermx) REAL ztsrf(ngridmx) REAL ztc(nlayermx), ztm(nlayermx+1) REAL zum(nlayermx+1) , zvm(nlayermx+1) REAL zqm(nlayermx+1,nqmx),zqm1(nlayermx+1) LOGICAL condsub(ngridmx) REAL subptimestep Integer Ntime real masse (nlayermx),w(nlayermx+1) real wq(ngridmx,nlayermx+1) real vstokes,reff real dWtotsn2 real condnconsn2(ngridmx) real nconsMAXn2 ! Special diagnostic variables real tconda1(ngridmx,nlayermx) real tconda2(ngridmx,nlayermx) real zdtsig (ngridmx,nlayermx) real zdtlatent (ngridmx,nlayermx) real zdt (ngridmx,nlayermx) REAL albediceF(ngridmx) SAVE albediceF INTEGER nsubtimestep,itsub !number of subtimestep when calling vl1d REAL subtimestep !ptimestep/nsubtimestep REAL dtmax REAL zplevhist(ngridmx) REAL zplevnew(ngridmx) REAL zplev(ngridmx) REAL zpicen2(ngridmx) REAL ztsrfhist(ngridmx) REAL zdtsrf(ngridmx) REAL globzplevnew REAL vmrn2(ngridmx) SAVE vmrn2 REAL stephan DATA stephan/5.67e-08/ ! Stephan Boltzman constant SAVE stephan !----------------------------------------------------------------------- ! Saved local variables ! REAL latcond REAL ccond REAL cpice ! for atm condensation SAVE cpice ! SAVE latcond,ccond SAVE ccond LOGICAL firstcall SAVE firstcall REAL SSUM EXTERNAL SSUM ! DATA latcond /2.5e5/ ! DATA latcond /1.98e5/ DATA cpice /1300./ DATA firstcall/.true./ INTEGER,SAVE :: i_n2ice=0 ! n2 ice CHARACTER(LEN=20) :: tracername ! to temporarily store text logical olkin ! option to prevent N2 ice effect in the south DATA olkin/.false./ save olkin CHARACTER(len=10) :: tname !----------------------------------------------------------------------- ! Initialisation IF (firstcall) THEN ccond=cpp/(g*lw_n2) print*,'In condens_n2cloud: ccond=',ccond,' latcond=',lw_n2 ! calculate global mean surface pressure for the fast mode IF (fast) THEN DO ig=1,ngridmx kp(ig) = exp(-phisfi(ig)/(r*38.)) ENDDO p00=glob_average2d(kp) ! mean pres at ref level ENDIF vmrn2(:) = 1. IF (ch4lag) then DO ig=1,ngridmx if (lati(ig)*180./pi.ge.latlag) then vmrn2(ig) = vmrlag endif ENDDO ENDIF firstcall=.false. ENDIF !----------------------------------------------------------------------- ! Calculation of n2 condensation / sublimation ! Variables used: ! picen2(ngrid) : amount of n2 ice on the ground (kg/m2) ! zcondicea(ngrid,nlayermx): condensation rate in layer l (kg/m2/s) ! zcondices(ngrid) : condensation rate on the ground (kg/m2/s) ! zfallice(ngrid,nlayermx) : amount of ice falling from layer l (kg/m2/s) ! zdtlatent(ngrid,nlayermx): dT/dt due to phase changes (K/s) ! Tendencies initially set to 0 zcondices(1:ngrid) = 0. pdtsrfc(1:ngrid) = 0. pdpsrf(1:ngrid) = 0. ztsrfhist(1:ngrid) = 0. condsub(1:ngrid) = .false. pdicen2(1:ngrid) = 0. zfallheat=0 pdqc(1:ngrid,1:nlayer,1:nq)=0 pdtc(1:ngrid,1:nlayer)=0 pduc(1:ngrid,1:nlayer)=0 pdvc(1:ngrid,1:nlayer)=0 zfallice(1:ngrid,1:nlayer+1)=0 zcondicea(1:ngrid,1:nlayer)=0 zdtlatent(1:ngrid,1:nlayer)=0 zt(1:ngrid,1:nlayer)=0. !----------------------------------------------------------------------- ! Atmospheric condensation ! Condensation / sublimation in the atmosphere ! -------------------------------------------- ! (calcul of zcondicea , zfallice and pdtc) zt(1:ngrid,1:nlayer)=pt(1:ngrid,1:nlayer)+ pdt(1:ngrid,1:nlayer)*ptimestep if (igcm_n2.ne.0) then zqn2(1:ngrid,1:nlayer) = 1. ! & temporaire ! zqn2(1:ngrid,1:nlayer)= pq(1:ngrid,1:nlayer,igcm_n2) + pdq(1:ngrid,1:nlayer,igcm_n2)*ptimestep else zqn2(1:ngrid,1:nlayer) = 1. end if if (.not.fast) then ! Forecast the atmospheric frost temperature 'ztcond' with function tcond_n2 DO l=1,nlayer DO ig=1,ngrid ztcond (ig,l) = tcond_n2(pplay(ig,l),zqn2(ig,l)) ENDDO ENDDO DO l=nlayer,1,-1 DO ig=1,ngrid pdtc(ig,l)=0. ! final tendancy temperature set to 0 IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN condsub(ig)=.true. !condensation at level l IF (zfallice(ig,l+1).gt.0) then zfallheat=zfallice(ig,l+1)*& (pphi(ig,l+1)-pphi(ig,l) +& cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/lw_n2 ELSE zfallheat=0. ENDIF zdtlatent(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))& *ccond*zdtlatent(ig,l)- zfallheat ! Case when the ice from above sublimes entirely ! """"""""""""""""""""""""""""""""""""""""""""""" IF ((zfallice(ig,l+1).lt.-zcondicea(ig,l)) & .AND. (zfallice(ig,l+1).gt.0)) THEN zdtlatent(ig,l)=(-zfallice(ig,l+1)+zfallheat)/& (ccond*(pplev(ig,l)-pplev(ig,l+1))) zcondicea(ig,l)= -zfallice(ig,l+1) END IF zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1) END IF ENDDO ENDDO endif ! not fast !----------------------------------------------------------------------- ! Condensation/sublimation on the ground ! (calculation of zcondices and pdtsrfc) ! Adding subtimesteps : in fast version, pressures too low lead to ! instabilities. IF (fast) THEN IF (pplev(1,1).gt.0.3) THEN nsubtimestep= 1 ELSE nsubtimestep= nbsub !max(nint(ptimestep/dtmax),1) ENDIF ELSE nsubtimestep= 1 ! if more then kp and p00 have to be calculated ! for nofast mode ENDIF subtimestep=ptimestep/float(nsubtimestep) do itsub=1,nsubtimestep ! first loop : getting zplev, ztsurf IF (itsub.eq.1) then DO ig=1,ngrid zplev(ig)=pplev(ig,1) ztsrfhist(ig)=ptsrf(ig) + pdtsrf(ig)*ptimestep ztsrf(ig)=ptsrf(ig) + pdtsrf(ig)*subtimestep !! zpicen2(ig)=picen2(ig) ENDDO ELSE ! additional loop : ! 1) pressure updated ! 2) direct redistribution of pressure over the globe ! 3) modification pressure for unstable cases ! 4) pressure update to conserve mass ! 5) temperature updated with radiative tendancies DO ig=1,ngrid zplevhist(ig)=zplev(ig) zplevnew(ig)=zplev(ig)+pdpsrf(ig)*subtimestep ! 1) !IF (zplevnew(ig).le.0.0001) then ! zplevnew(ig)=0.0001*kp(ig)/p00 !ENDIF ENDDO ! intermediaire de calcul: valeur moyenne de zplevnew (called twice in the code) globzplevnew=glob_average2d(zplevnew) DO ig=1,ngrid zplev(ig)=kp(ig)*globzplevnew/p00 ! 2) ENDDO DO ig=1,ngrid ! 3) unstable case condensation and sublimation: zplev=zplevhist IF (((pdpsrf(ig).lt.0.).and.(tcond_n2(zplev(ig),zqn2(ig,1)).le.ztsrf(ig))).or. & ((pdpsrf(ig).gt.0.).and.(tcond_n2(zplev(ig),zqn2(ig,1)).ge.ztsrf(ig)))) then zplev(ig)=zplevhist(ig) ENDIF zplevhist(ig)=zplev(ig) ENDDO zplev=zplev*globzplevnew/glob_average2d(zplevhist) ! 4) DO ig=1,ngrid ! 5) zdtsrf(ig)=pdtsrf(ig) + (stephan/pcapcal(ig))*(ptsrf(ig)**4-ztsrf(ig)**4) ztsrf(ig)=ztsrf(ig)+pdtsrfc(ig)*subtimestep+zdtsrf(ig)*subtimestep zpicen2(ig)=zpicen2(ig)+pdicen2(ig)*subtimestep ENDDO ENDIF ! (itsub=1) DO ig=1,ngrid ! forecast of frost temperature ztcondsol !ztcondsol(ig) = tcond_n2(zplev(ig),zqn2(ig,1)) ztcondsol(ig) = tcond_n2(zplev(ig),vmrn2(ig)) ! Loop over where we have condensation / sublimation IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR. & ! ground cond ((ztsrf(ig) .GT. ztcondsol(ig)) .AND. & ! ground sublim (zpicen2(ig) .GT. 0.))) THEN condsub(ig) = .true. ! condensation or sublimation ! Condensation or partial sublimation of N2 ice if (ztsrf(ig) .LT. ztcondsol(ig)) then ! condensation ! Include a correction to account for the cooling of air near ! the surface before condensing: zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) & /((lw_n2+cpp*(zt(ig,1)-ztcondsol(ig)))*subtimestep) else ! sublimation zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) & /(lw_n2*subtimestep) end if pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/subtimestep ! partial sublimation of N2 ice ! If the entire N_2 ice layer sublimes ! (including what has just condensed in the atmosphere) IF((zpicen2(ig)/subtimestep).LE. & -zcondices(ig))THEN zcondices(ig) = -zpicen2(ig)/subtimestep pdtsrfc(ig)=(lw_n2/pcapcal(ig))* & (zcondices(ig)) END IF ! Changing N2 ice amount and pressure pdicen2(ig) = zcondices(ig) pdpsrf(ig) = -pdicen2(ig)*g ! pdpsrf(ig) = 0. ! OPTION to check impact N2 sub/cond IF (zplev(ig)+pdpsrf(ig)*subtimestep.le.0.0000001) then pdpsrf(ig)=(0.0000001*kp(ig)/p00-zplev(ig))/subtimestep pdicen2(ig)=-pdpsrf(ig)/g ENDIF ELSE ! no condsub pdpsrf(ig)=0. pdicen2(ig)=0. pdtsrfc(ig)=0. ENDIF ENDDO ! ig enddo ! subtimestep ! Updating pressure, temperature and ice reservoir DO ig=1,ngrid pdpsrf(ig)=(zplev(ig)+pdpsrf(ig)*subtimestep-pplev(ig,1))/ptimestep ! Two options here : 1 ok, 2 is wrong pdicen2(ig)=(zpicen2(ig)+pdicen2(ig)*subtimestep-picen2(ig))/ptimestep !pdicen2(ig)=-pdpsrf(ig)/g pdtsrfc(ig)=((ztsrf(ig)+pdtsrfc(ig)*subtimestep)-(ztsrfhist(ig)))/ptimestep ! security if (picen2(ig) + pdicen2(ig)*ptimestep.lt.0.) then write(*,*) 'WARNING in condense_n2:' write(*,*) picen2(ig),pdicen2(ig)*ptimestep pdicen2(ig)= -picen2(ig)/ptimestep pdpsrf(ig)=-pdicen2(ig)*g endif if(.not.picen2(ig).ge.0.) THEN ! if(picen2(ig) + pdicen2(ig)*ptimestep.le.-1.e-8) then print*, 'WARNING NEG RESERVOIR in condense_n2: picen2(',ig,')=', picen2(ig) + pdicen2(ig)*ptimestep ! pdicen2(ig)= -picen2(ig)/ptimestep ! else picen2(ig)=0.0 ! endif endif ENDDO ! *************************************************************** ! Correction to account for redistribution between sigma or hybrid ! layers when changing surface pressure (and warming/cooling ! of the n2 currently changing phase). ! ************************************************************* if (.not.fast) then DO ig=1,ngrid if (condsub(ig)) then ! Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up) ! """""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" zmflux(1) = -zcondices(ig) DO l=1,nlayer zmflux(l+1) = zmflux(l) -zcondicea(ig,l) & + (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1)) ! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0. END DO ! Mass of each layer ! ------------------ DO l=1,nlayer masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g END DO ! Corresponding fluxes for T,U,V,Q ! """""""""""""""""""""""""""""""" ! averaging operator for TRANSPORT ! """""""""""""""""""""""""""""""" ! Subtimestep loop to perform the redistribution gently and simultaneously with ! the other tendencies ! Estimation of subtimestep (using only the first layer, the most critical) dtmax=ptimestep if (zmflux(1).gt.1.e-20) then dtmax=min(dtmax,masse(1)*zqn2(ig,1)/abs(zmflux(1))) endif nsubtimestep= max(nint(ptimestep/dtmax),nint(2.)) subtimestep=ptimestep/float(nsubtimestep) ! New flux for each subtimestep do l=1,nlayer+1 w(l)=-zmflux(l)*subtimestep enddo ! initializing variables that will vary during subtimestep: do l=1,nlayer ztc(l) =pt(ig,l) zu(l) =pu(ig,l) zv(l) =pv(ig,l) do iq=1,nqmx zq(l,iq) = pq(ig,l,iq) enddo end do ! loop over nsubtimestep ! """""""""""""""""""""" do itsub=1,nsubtimestep ! Progressively adding tendancies from other processes. do l=1,nlayer ztc(l) =ztc(l) +(pdt(ig,l) + zdtlatent(ig,l))*subtimestep zu(l) =zu(l) +pdu( ig,l) * subtimestep zv(l) =zv(l) +pdv( ig,l) * subtimestep do iq=1,nqmx zq(l,iq) = zq(l,iq) + pdq(ig,l,iq)* subtimestep enddo end do ! Change of mass in each layer do l=1,nlayer masse(l)=masse(l)+pdpsrf(ig)*subtimestep*(pplev(ig,l) - pplev(ig,l+1))& /(g*pplev(ig,1)) end do ztm(2:nlayermx+1)=0. zum(2:nlayermx+1)=0. zvm(2:nlayermx+1)=0. zqm1(1:nlayermx+1)=0. ! Van Leer scheme: call vl1d(ztc,2.,masse,w,ztm) call vl1d(zu ,2.,masse,w,zum) call vl1d(zv ,2.,masse,w,zvm) do iq=1,nqmx do l=1,nlayer zq1(l)=zq(l,iq) enddo zqm1(1)=zqm(1,iq) call vl1d(zq1,2.,masse,w,zqm1) do l=2,nlayer zqm(l,iq)=zqm1(l) enddo enddo ! Correction to prevent negative value for qn2 if (igcm_n2.ne.0) then do l=1,nlayer-1 if (w(l)*zqm(l,igcm_n2).gt.zq(l,igcm_n2)*masse(l)) then zqm(l+1,igcm_n2)=max(zqm(l+1,igcm_n2), & (zqm(l,igcm_n2)*w(l) -zq(l,igcm_n2)*masse(l))/w(l+1) ) else exit endif end do end if ! Value transfert at the surface interface when condensation sublimation: if (zmflux(1).lt.0) then ! Surface condensation zum(1)= zu(1) zvm(1)= zv(1) ztm(1) = ztc(1) else ! Surface sublimation: ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep zum(1) = 0 zvm(1) = 0 end if do iq=1,nqmx zqm(1,iq)=0. ! most tracer do not condense ! enddo ! Special case if the tracer is n2 gas if (igcm_n2.ne.0) zqm(1,igcm_n2)=1. !!! Source haze: 0.02 pourcent when n2 sublimes IF (source_haze) THEN IF (pdicen2(ig).lt.0) THEN DO iq=1,nq tname=noms(iq) if (tname(1:4).eq."haze") then !zqm(1,iq)=0.02 !zqm(1,iq)=-pdicen2(ig)*0.02 zqm(1,iq)=-pdicen2(ig)*ptimestep*0.02 !zqm(10,iq)=-pdicen2(ig)*ptimestep*100. !zqm(1,iq)=-pdicen2(ig)*1000000. endif ENDDO ENDIF ENDIF ztm(nlayer+1)= ztc(nlayer) ! should not be used, but... zum(nlayer+1)= zu(nlayer) ! should not be used, but... zvm(nlayer+1)= zv(nlayer) ! should not be used, but... do iq=1,nqmx zqm(nlayer+1,iq)= zq(nlayer,iq) enddo ! Tendencies on T, U, V, Q ! """"""""""""""""""""""" DO l=1,nlayer ! Tendencies on T zdtsig(ig,l) = (1/masse(l)) * & ( zmflux(l)*(ztm(l) - ztc(l)) & - zmflux(l+1)*(ztm(l+1) - ztc(l)) & + zcondicea(ig,l)*(ztcond(ig,l)-ztc(l)) ) ! Tendencies on U pduc(ig,l) = (1/masse(l)) * & ( zmflux(l)*(zum(l) - zu(l))& - zmflux(l+1)*(zum(l+1) - zu(l)) ) ! Tendencies on V pdvc(ig,l) = (1/masse(l)) * & ( zmflux(l)*(zvm(l) - zv(l)) & - zmflux(l+1)*(zvm(l+1) - zv(l)) ) END DO ! Tendencies on Q do iq=1,nqmx if (iq.eq.igcm_n2) then ! SPECIAL Case when the tracer IS N2 : DO l=1,nlayer pdqc(ig,l,iq)= (1/masse(l)) * & ( zmflux(l)*(zqm(l,iq) - zq(l,iq)) & - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))& + zcondicea(ig,l)*(zq(l,iq)-1.) ) END DO else DO l=1,nlayer pdqc(ig,l,iq)= (1/masse(l)) * & ( zmflux(l)*(zqm(l,iq) - zq(l,iq)) & - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq)) & + zcondicea(ig,l)*zq(l,iq) ) END DO end if enddo ! Update variables at the end of each subtimestep. do l=1,nlayer ztc(l) =ztc(l) + zdtsig(ig,l) *subtimestep zu(l) =zu(l) + pduc(ig,l) *subtimestep zv(l) =zv(l) + pdvc(ig,l) *subtimestep do iq=1,nqmx zq(l,iq) = zq(l,iq) + pdqc(ig,l,iq)*subtimestep enddo end do enddo ! loop on nsubtimestep ! Recomputing Total tendencies do l=1,nlayer pdtc(ig,l) = (ztc(l) - zt(ig,l) )/ptimestep pduc(ig,l) = (zu(l) - (pu(ig,l) + pdu(ig,l)*ptimestep))/ptimestep pdvc(ig,l) = (zv(l) - (pv(ig,l) + pdv(ig,l)*ptimestep))/ptimestep do iq=1,nqmx pdqc(ig,l,iq) = (zq(l,iq) - (pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep))/ptimestep ! Correction temporaire if (iq.eq.igcm_n2) then if((pq(ig,l,iq) +(pdqc(ig,l,iq)+ pdq(ig,l,iq))*ptimestep) & .lt.0.01) then ! if n2 < 1 % ! pdqc(ig,l,iq)=(0.01-pq(ig,l,iq))/ptimestep-pdq(ig,l,iq) end if end if enddo end do ! *******************************TEMPORAIRE ****************** if (ngrid.eq.1) then write(*,*) 'nsubtimestep=' ,nsubtimestep write(*,*) 'masse avant' , (pplev(ig,1) - pplev(ig,2))/g write(*,*) 'masse apres' , masse(1) write(*,*) 'zmflux*DT, l=1' , zmflux(1)*ptimestep write(*,*) 'zmflux*DT, l=2' , zmflux(2)*ptimestep write(*,*) 'pq, l=1,2,3' , pq(1,1,1), pq(1,2,1),pq(1,3,1) write(*,*) 'zq, l=1,2,3' , zq(1,1), zq(2,1),zq(3,1) write(*,*) 'dq*Dt l=1' , pdq(1,1,1)*ptimestep write(*,*) 'dqcond*Dt l=1' , pdqc(1,1,1)*ptimestep end if ! *********************************************************** end if ! if (condsub) END DO ! loop on ig endif ! not fast ! ************ Option Olkin to prevent N2 effect in the south******** 112 continue if (olkin) then DO ig=1,ngrid if (lati(ig).lt.0) then pdtsrfc(ig) = max(0.,pdtsrfc(ig)) pdpsrf(ig) = 0. pdicen2(ig) = 0. do l=1,nlayer pdtc(ig,l) = max(0.,zdtlatent(ig,l)) pduc(ig,l) = 0. pdvc(ig,l) = 0. do iq=1,nqmx pdqc(ig,l,iq) = 0 enddo end do end if END DO end if ! ******************************************************************* ! *************************************************************** ! Ecriture des diagnostiques ! *************************************************************** ! DO l=1,nlayer ! DO ig=1,ngrid ! Taux de cond en kg.m-2.pa-1.s-1 ! tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1)) ! Taux de cond en kg.m-3.s-1 ! tconda2(ig,l)=tconda1(ig,l)*pplay(ig,l)*g/(r*pt(ig,l)) ! END DO ! END DO ! call WRITEDIAGFI(ngridmx,'tconda1', & ! 'Taux de condensation N2 atmospherique /Pa', & ! 'kg.m-2.Pa-1.s-1',3,tconda1) ! call WRITEDIAGFI(ngridmx,'tconda2', & ! 'Taux de condensation N2 atmospherique /m', & ! 'kg.m-3.s-1',3,tconda2) return end subroutine condens_n2 !------------------------------------------------------------------------- real function tcond_n2(p,vmr) ! Calculates the condensation temperature for N2 at pressure P and vmr implicit none real, intent(in):: p,vmr ! tcond_n2 = (1.)/(0.026315-0.0011877*log(.7143*p*vmr)) IF (p.ge.0.529995) then ! tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB ! tcond_n2 = (1.)/(1./63.1470-296.925/(2.5e5*0.98)*log(1./(0.125570*1.e5)*p*vmr)) tcond_n2 = (1.)/(0.01583606505-1.211938776e-3*log(7.963685594e-5*p*vmr)) ELSE ! tcond Fray and Schmitt for N2 phase alpha(T<35.6 K) FIT BT ! tcond_n2 = (1.)/(1./35.6-296.925/(2.5e5*1.09)*log(1./(0.508059)*p*vmr)) tcond_n2 = (1.)/(1./35.6-1.089633028e-3*log(1.968275338*p*vmr)) ENDIF return end function tcond_n2 !------------------------------------------------------------------------- real function pcond_n2(t,vmr) ! Calculates the condensation pressure for N2 at temperature T and vmr implicit none real, intent(in):: t,vmr ! tcond_n2 = (1.)/(0.026315-0.0011877*log(.7143*p*vmr)) IF (t.ge.35.6) then ! tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB ! pcond_n2 = 0.125570*1.e5/vmr*exp((2.5e5*0.98)/296.925*(1./63.1470-1./t)) pcond_n2 = 0.125570e5/vmr*exp(825.1241896*(1./63.147-1./t)) ELSE ! tcond Fray and Schmitt for N2 phase alpha(T<35.6 K) FIT TB ! pcond_n2 = 0.508059/vmr*exp((2.5e5*1.09)/296.925*(1./35.6-1./t)) pcond_n2 = 0.508059/vmr*exp(917.7401701*(1./35.6-1./t)) ENDIF return end function pcond_n2 !------------------------------------------------------------------------- real function glob_average2d(var) ! Calculates the global average of variable var use comgeomfi_h implicit none #include "dimensions.h" #include "dimphys.h" ! INTEGER ngrid REAL var(ngridmx) INTEGER ig glob_average2d = 0. DO ig=2,ngridmx-1 glob_average2d = glob_average2d + var(ig)*area(ig) END DO glob_average2d = glob_average2d + var(1)*area(1)*iim glob_average2d = glob_average2d + var(ngridmx)*area(ngridmx)*iim glob_average2d = glob_average2d/(totarea+(area(1)+area(ngridmx))*(iim-1)) end function glob_average2d ! ***************************************************************** subroutine vl1d(q,pente_max,masse,w,qm) ! ! Operateur de moyenne inter-couche pour calcul de transport type ! Van-Leer " pseudo amont " dans la verticale ! q,w sont des arguments d'entree pour le s-pg .... ! masse : masse de la couche Dp/g ! w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2) ! pente_max = 2 conseillee ! -------------------------------------------------------------------- IMPLICIT NONE #include "dimensions.h" ! Arguments: ! ---------- real masse(llm),pente_max REAL q(llm),qm(llm+1) REAL w(llm+1) ! ! Local ! --------- ! INTEGER l ! real dzq(llm),dzqw(llm),adzqw(llm),dzqmax real sigw, Mtot, MQtot integer m ! On oriente tout dans le sens de la pression ! W > 0 WHEN DOWN !!!!!!!!!!!!! do l=2,llm dzqw(l)=q(l-1)-q(l) adzqw(l)=abs(dzqw(l)) enddo do l=2,llm-1 if(dzqw(l)*dzqw(l+1).gt.0.) then dzq(l)=0.5*(dzqw(l)+dzqw(l+1)) else dzq(l)=0. endif dzqmax=pente_max*min(adzqw(l),adzqw(l+1)) dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l)) enddo dzq(1)=0. dzq(llm)=0. do l = 1,llm-1 ! Regular scheme (transfered mass < layer mass) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then sigw=w(l+1)/masse(l+1) qm(l+1)=(q(l+1)+0.5*(1.-sigw)*dzq(l+1)) else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then sigw=w(l+1)/masse(l) qm(l+1)=(q(l)-0.5*(1.+sigw)*dzq(l)) ! Extended scheme (transfered mass > layer mass) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ else if(w(l+1).gt.0.) then m=l+1 Mtot = masse(m) MQtot = masse(m)*q(m) do while ((m.lt.llm).and.(w(l+1).gt.(Mtot+masse(m+1)))) m=m+1 Mtot = Mtot + masse(m) MQtot = MQtot + masse(m)*q(m) end do if (m.lt.llm) then sigw=(w(l+1)-Mtot)/masse(m+1) qm(l+1)= (1/w(l+1))*(MQtot + (w(l+1)-Mtot)* & (q(m+1)+0.5*(1.-sigw)*dzq(m+1)) ) else w(l+1) = Mtot qm(l+1) = Mqtot / Mtot write(*,*) 'top layer is disapearing !' stop end if else ! if(w(l+1).lt.0) m = l-1 Mtot = masse(m+1) MQtot = masse(m+1)*q(m+1) do while((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(max(m,1))))) m=m-1 Mtot = Mtot + masse(m+1) MQtot = MQtot + masse(m+1)*q(m+1) end do if (m.gt.0) then sigw=(w(l+1)+Mtot)/masse(m) qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)* & (q(m)-0.5*(1.+sigw)*dzq(m)) ) else qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)*qm(1)) end if end if enddo return end subroutine vl1d