subroutine condense_n2(klon,klev,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 use comcstfi_mod, only: g, r, cpp, pi USE surfdat_h, only: phisfi,kp,p00 USE tracer_h, only: noms, igcm_n2, lw_n2 USE callkeys_mod, only: fast,ch4lag,latlag,nbsub,no_n2frost,tsurfmax,kmixmin,source_haze,vmrlag USE comvert_mod, ONLY: ap,bp use geometry_mod, only: latitude implicit none !================================================================== ! Purpose ! ------- ! Condense and/or sublime N2 ice on the ground and in the ! atmosphere, and sediment the ice. ! ! Inputs ! ------ ! klon Number of vertical columns ! klev Number of layers ! pplay(klon,klev) Pressure layers ! pplev(klon,klev+1) Pressure levels ! pt(klon,klev) Temperature (in K) ! ptsrf(klon) Surface temperature ! ! pdt(klon,klev) Time derivative before condensation/sublimation of pt ! pdtsrf(klon) Time derivative before condensation/sublimation of ptsrf ! picen2(klon) n2 ice at the surface (kg/m2) ! ! Outputs ! ------- ! pdpsrf(klon) \ Contribution of condensation/sublimation ! pdtc(klon,klev) / to the time derivatives of Ps, pt, and ptsrf ! pdtsrfc(klon) / ! pdicen2(klon) Tendancy n2 ice at the surface (kg/m2) ! ! Both ! ---- ! psolaralb(klon) Albedo at the surface ! pemisurf(klon) Emissivity of the surface ! ! Authors ! ------- ! Francois Forget (1996,2013) ! Converted to Fortran 90 and slightly modified by R. Wordsworth (2009) ! ! !================================================================== !----------------------------------------------------------------------- ! Arguments INTEGER klon, klev, nq REAL ptimestep REAL pplay(klon,klev),pplev(klon,klev+1) REAL pcapcal(klon) REAL pt(klon,klev) REAL ptsrf(klon),flu1(klon),flu2(klon),flu3(klon) REAL pphi(klon,klev) REAL pdt(klon,klev),pdtsrf(klon),pdtc(klon,klev) REAL pdtsrfc(klon),pdpsrf(klon) REAL picen2(klon),psolaralb(klon),pemisurf(klon) REAL pu(klon,klev) , pv(klon,klev) REAL pdu(klon,klev) , pdv(klon,klev) REAL pduc(klon,klev) , pdvc(klon,klev) REAL pq(klon,klev,nq),pdq(klon,klev,nq) REAL pdqc(klon,klev,nq) !----------------------------------------------------------------------- ! Local variables INTEGER l,ig,ilay,it,iq,ich4_gas REAL*8 zt(klon,klev) REAL tcond_n2 REAL pcond_n2 REAL glob_average2d ! function REAL zqn2(klon,klev) ! N2 MMR used to compute Tcond/zqn2 REAL ztcond (klon,klev) REAL ztcondsol(klon),zfallheat REAL pdicen2(klon) REAL zcondicea(klon,klev), zcondices(klon) REAL zfallice(klon,klev+1) REAL zmflux(klev+1) REAL zu(klev),zv(klev) REAL zq(klev,nq),zq1(klev) REAL ztsrf(klon) REAL ztc(klev), ztm(klev+1) REAL zum(klev+1) , zvm(klev+1) REAL zqm(klev+1,nq),zqm1(klev+1) LOGICAL condsub(klon) REAL subptimestep Integer Ntime real masse (klev),w(klev+1) real wq(klon,klev+1) real vstokes,reff real dWtotsn2 real condnconsn2(klon) real nconsMAXn2 ! Special diagnostic variables real tconda1(klon,klev) real tconda2(klon,klev) real zdtsig (klon,klev) real zdtlatent (klon,klev) real zdt (klon,klev) ! REAL albediceF(klon) ! SAVE albediceF INTEGER nsubtimestep,itsub !number of subtimestep when calling vl1d REAL subtimestep !ptimestep/nsubtimestep REAL dtmax REAL zplevhist(klon) REAL zplevnew(klon) REAL zplev(klon) REAL zpicen2(klon) REAL ztsrfhist(klon) REAL zdtsrf(klon) REAL globzplevnew real,dimension(:),save,allocatable :: vmrn2 !$OMP THREADPRIVATE(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 condense_n2cloud: ccond=',ccond,' latcond=',lw_n2 ! calculate global mean surface pressure for the fast mode IF (.not. ALLOCATED(kp)) ALLOCATE(kp(klon)) DO ig=1,klon kp(ig) = exp(-phisfi(ig)/(r*38.)) ENDDO IF (fast) THEN p00=glob_average2d(kp) ! mean pres at ref level ENDIF ALLOCATE(vmrn2(klon)) vmrn2(:) = 1. !IF (ch4lag) then ! DO ig=1,klon ! if (latitude(ig)*180./pi.ge.latlag) then ! vmrn2(ig) = vmrlag ! endif ! ENDDO !ENDIF IF (no_n2frost) then DO ig=1,klon if (picen2(ig).eq.0.) then vmrn2(ig) = 1.e-15 endif ENDDO ENDIF firstcall=.false. ENDIF !----------------------------------------------------------------------- ! Calculation of n2 condensation / sublimation ! Variables used: ! picen2(klon) : amount of n2 ice on the ground (kg/m2) ! zcondicea(klon,klev): condensation rate in layer l (kg/m2/s) ! zcondices(klon) : condensation rate on the ground (kg/m2/s) ! zfallice(klon,klev) : amount of ice falling from layer l (kg/m2/s) ! zdtlatent(klon,klev): dT/dt due to phase changes (K/s) ! Tendencies initially set to 0 zcondices(1:klon) = 0. pdtsrfc(1:klon) = 0. pdpsrf(1:klon) = 0. ztsrfhist(1:klon) = 0. condsub(1:klon) = .false. pdicen2(1:klon) = 0. zfallheat=0 pdqc(1:klon,1:klev,1:nq)=0 pdtc(1:klon,1:klev)=0 pduc(1:klon,1:klev)=0 pdvc(1:klon,1:klev)=0 zfallice(1:klon,1:klev+1)=0 zcondicea(1:klon,1:klev)=0 zdtlatent(1:klon,1:klev)=0 zt(1:klon,1:klev)=0. !----------------------------------------------------------------------- ! Atmospheric condensation ! Condensation / sublimation in the atmosphere ! -------------------------------------------- ! (calcul of zcondicea , zfallice and pdtc) zt(1:klon,1:klev)=pt(1:klon,1:klev)+ pdt(1:klon,1:klev)*ptimestep if (igcm_n2.ne.0) then zqn2(1:klon,1:klev) = 1. ! & temporaire ! zqn2(1:klon,1:klev)= pq(1:klon,1:klev,igcm_n2) + pdq(1:klon,1:klev,igcm_n2)*ptimestep else zqn2(1:klon,1:klev) = 1. end if if (.not.fast) then ! Forecast the atmospheric frost temperature 'ztcond' with function tcond_n2 DO l=1,klev DO ig=1,klon ztcond (ig,l) = tcond_n2(pplay(ig,l),zqn2(ig,l)) ENDDO ENDDO DO l=klev,1,-1 DO ig=1,klon 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,klon 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,klon 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,klon zplev(ig)=kp(ig)*globzplevnew/p00 ! 2) ENDDO DO ig=1,klon ! 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,klon ! 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,klon ! 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 (fast.and.(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,klon 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,klon 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,klev 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,klev 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,klev+1 w(l)=-zmflux(l)*subtimestep enddo ! initializing variables that will vary during subtimestep: do l=1,klev ztc(l) =pt(ig,l) zu(l) =pu(ig,l) zv(l) =pv(ig,l) do iq=1,nq 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,klev 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,nq zq(l,iq) = zq(l,iq) + pdq(ig,l,iq)* subtimestep enddo end do ! Change of mass in each layer do l=1,klev masse(l)=masse(l)+pdpsrf(ig)*subtimestep*(pplev(ig,l) - pplev(ig,l+1))& /(g*pplev(ig,1)) end do ! 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,nq 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. ztm(2:klev+1)=0. zum(2:klev+1)=0. zvm(2:klev+1)=0. zqm1(1:klev+1)=0. ! Van Leer scheme: call vl1d(klev,ztc,2.,masse,w,ztm) call vl1d(klev,zu ,2.,masse,w,zum) call vl1d(klev,zv ,2.,masse,w,zvm) do iq=1,nq do l=1,klev zq1(l)=zq(l,iq) enddo zqm1(1)=zqm(1,iq) call vl1d(klev,zq1,2.,masse,w,zqm1) do l=2,klev zqm(l,iq)=zqm1(l) enddo enddo ! Correction to prevent negative value for qn2 if (igcm_n2.ne.0) then zqm(1,igcm_n2)=1. do l=1,klev-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,nq 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(klev+1)= ztc(klev) ! should not be used, but... zum(klev+1)= zu(klev) ! should not be used, but... zvm(klev+1)= zv(klev) ! should not be used, but... do iq=1,nq zqm(klev+1,iq)= zq(klev,iq) enddo ! Tendencies on T, U, V, Q ! """"""""""""""""""""""" DO l=1,klev ! 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,nq if (iq.eq.igcm_n2) then ! SPECIAL Case when the tracer IS N2 : DO l=1,klev 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,klev 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,klev 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,nq zq(l,iq) = zq(l,iq) + pdqc(ig,l,iq)*subtimestep enddo end do enddo ! loop on nsubtimestep ! Recomputing Total tendencies do l=1,klev 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,nq 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 % ! write(*,*) 'Warning: 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 (klon.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,klon if (latitude(ig).lt.0) then pdtsrfc(ig) = max(0.,pdtsrfc(ig)) pdpsrf(ig) = 0. pdicen2(ig) = 0. do l=1,klev pdtc(ig,l) = max(0.,zdtlatent(ig,l)) pduc(ig,l) = 0. pdvc(ig,l) = 0. do iq=1,nq pdqc(ig,l,iq) = 0 enddo end do end if END DO end if ! ******************************************************************* ! *************************************************************** ! Ecriture des diagnostiques ! *************************************************************** ! DO l=1,klev ! DO ig=1,klon ! 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(klon,'tconda1', & ! 'Taux de condensation N2 atmospherique /Pa', & ! 'kg.m-2.Pa-1.s-1',3,tconda1) ! call WRITEDIAGFI(klon,'tconda2', & ! 'Taux de condensation N2 atmospherique /m', & ! 'kg.m-3.s-1',3,tconda2) return end subroutine condense_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 use dimphy, only: klon USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat use geometry_mod, only: cell_area, latitude implicit none ! INTEGER klon REAL var(klon) INTEGER ig glob_average2d = 0. DO ig=2,klon-1 glob_average2d = glob_average2d + var(ig)*cell_area(ig) END DO glob_average2d = glob_average2d + var(1)*cell_area(1)*nbp_lon glob_average2d = glob_average2d + var(klon)*cell_area(klon)*nbp_lon glob_average2d = glob_average2d/(totarea+(cell_area(1)+cell_area(klon))*(nbp_lon-1)) end function glob_average2d ! ***************************************************************** subroutine vl1d(klev,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 ! Arguments: ! ---------- integer klev real masse(klev),pente_max REAL q(klev),qm(klev+1) REAL w(klev+1) ! ! Local ! --------- ! INTEGER l ! real dzq(klev),dzqw(klev),adzqw(klev),dzqmax real sigw, Mtot, MQtot integer m ! On oriente tout dans le sens de la pression ! W > 0 WHEN DOWN !!!!!!!!!!!!! do l=2,klev dzqw(l)=q(l-1)-q(l) adzqw(l)=abs(dzqw(l)) enddo do l=2,klev-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(klev)=0. do l = 1,klev-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) if (m.lt.klev) then ! because some compilers will have problems ! evaluating masse(klev+1) do while ((m.lt.klev).and.(w(l+1).gt.(Mtot+masse(m+1)))) m=m+1 Mtot = Mtot + masse(m) MQtot = MQtot + masse(m)*q(m) if (m.eq.klev) exit end do endif if (m.lt.klev) 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) if (m.gt.0) then ! because some compilers will have problems ! evaluating masse(0) do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m)))) m=m-1 Mtot = Mtot + masse(m+1) MQtot = MQtot + masse(m+1)*q(m+1) if (m.eq.0) exit end do endif 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