subroutine moistadj_generic(ngrid, nlayer, nq, pt, pq, pdq, pplev, pplay, pdtmana, pdqmana, ptimestep, rneb) use generic_cloud_common_h use generic_tracer_index_mod, only: generic_tracer_index use tracer_h use ioipsl_getin_p_mod, only: getin_p !-> to get the metallicity use comcstfi_mod, only: r, cpp, mugaz use callkeys_mod, only: moist_convection_inhibition implicit none !===================================================================== ! ! Purpose ! ------- ! Calculates moist convective adjustment by the method of Manabe. ! ! Authors ! ------- ! Adapted from the moistadj.F90 routine ! for generic condensable species (GCS) tracers ! by Noe CLEMENT (2023) ! !===================================================================== INTEGER,INTENT(IN) :: ngrid, nlayer, nq REAL,INTENT(IN) :: pt(ngrid,nlayer) ! temperature (K) REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracer (kg/kg) REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq) REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa) REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa) REAL,INTENT(IN) :: ptimestep ! physics timestep (s) REAL,INTENT(OUT) :: pdqmana(ngrid,nlayer,nq) ! tracer tendencies (kg/kg.s-1) REAL,INTENT(OUT) :: pdtmana(ngrid,nlayer) ! temperature increment(K/s) REAL,INTENT(OUT) :: rneb(ngrid,nlayer) ! cloud fraction ! Options : !real, save :: metallicity ! metallicity of planet REAL,SAVE :: metallicity = 0.0 !$OMP THREADPRIVATE(metallicity) ! local variables REAL zt(ngrid,nlayer) ! temperature (K) REAL zq(ngrid,nlayer) ! humidite specifique (kg/kg) REAL d_t(ngrid,nlayer) ! temperature increment REAL d_q(ngrid,nlayer) ! incrementation pour vapeur d'eau REAL d_ql(ngrid,nlayer) ! incrementation pour l'eau liquide ! REAL t_coup ! PARAMETER (t_coup=234.0) REAL seuil_vap PARAMETER (seuil_vap=1.0E-10) ! Local variables logical call_ice_vap_generic ! to call only one time the ice/vap pair of a tracer integer igcm_generic_vap, igcm_generic_ice ! index of the vap and ice of generic_tracer INTEGER i, k, iq, nn INTEGER, PARAMETER :: niter = 1 INTEGER k1, k1p, k2, k2p LOGICAL itest(ngrid) REAL delta_q(ngrid, nlayer) DOUBLE PRECISION :: cp_new_t(nlayer), v_cptt(ngrid,nlayer) REAL cp_delta_t(nlayer) DOUBLE PRECISION :: v_cptj(nlayer), v_cptjk1, v_ssig REAL v_p, v_t, v_zqs,v_cptt2,v_pratio,v_dlnpsat REAL zqs(ngrid,nlayer), zdqs(ngrid,nlayer),zpsat(ngrid,nlayer),zdlnpsat(ngrid,nlayer) REAL zq1(ngrid), zq2(ngrid) DOUBLE PRECISION :: gamcpdz(ngrid,2:nlayer) DOUBLE PRECISION :: zdp, zdpm real q_cri(ngrid,nlayer) ! moist convection inhibition criterion REAL zsat ! super-saturation REAL zflo ! flotabilite DOUBLE PRECISION :: local_q(ngrid,nlayer),local_t(ngrid,nlayer) REAL zdelta, zcor, zcvm5 REAL dEtot, dqtot, masse ! conservation diagnostics real dL1tot, dL2tot ! Indices of generic vapour and ice tracers real,save :: RCPD=0.0 INTEGER,SAVE :: i_vap_generic=0 ! Generic Condensable Species vapour INTEGER,SAVE :: i_ice_generic=0 ! Generic Condensable Species ice !$OMP THREADPRIVATE(i_vap_generic,i_ice_generic,RCPD) LOGICAL,SAVE :: firstcall=.TRUE. !$OMP THREADPRIVATE(firstcall) IF (firstcall) THEN RCPD = cpp write(*,*) "value for metallicity? " metallicity=0.0 ! default value call getin_p("metallicity",metallicity) write(*,*) " metallicity = ",metallicity 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 ! to call only one time the ice/vap pair of a tracer write(*,*) "moistadjustment : taking back the values you have set in 'table_tracers_condensable' for thermodynamics. If you have forgotten one, it will crash." m = constants_mass(iq) delta_vapH = constants_delta_vapH(iq) epsi_generic = constants_epsi_generic(iq) RLVTT_generic = constants_RLVTT_generic(iq) RCPV_generic = constants_RCPV_generic(iq) Tref = constants_Tref(iq) write(*,*) noms(igcm_generic_vap),", q_cri at ", Tref, "K (in kg/kg): ", ( 1 / (1 - 1/epsi_generic)) * (r * mugaz/1000.) / delta_vapH * Tref endif enddo i_vap_generic=igcm_generic_vap i_ice_generic=igcm_generic_ice write(*,*) "rain: i_ice_generic=",i_ice_generic write(*,*) " i_vap_generic=",i_vap_generic firstcall = .FALSE. ENDIF ! GCM -----> subroutine variables zq(1:ngrid,1:nlayer) = pq(1:ngrid,1:nlayer,i_vap_generic)+ pdq(1:ngrid,1:nlayer,i_vap_generic)*ptimestep zt(1:ngrid,1:nlayer) = pt(1:ngrid,1:nlayer) pdqmana(1:ngrid,1:nlayer,1:nq)=0.0 DO k = 1, nlayer DO i = 1, ngrid if(zq(i,k).lt.0.)then zq(i,k)=0.0 endif ENDDO ENDDO local_q(1:ngrid,1:nlayer) = zq(1:ngrid,1:nlayer) local_t(1:ngrid,1:nlayer) = zt(1:ngrid,1:nlayer) rneb(1:ngrid,1:nlayer) = 0.0 d_ql(1:ngrid,1:nlayer) = 0.0 d_t(1:ngrid,1:nlayer) = 0.0 d_q(1:ngrid,1:nlayer) = 0.0 ! Calculate v_cptt DO k = 1, nlayer DO i = 1, ngrid v_cptt(i,k) = RCPD * local_t(i,k) v_t = MAX(local_t(i,k),15.) v_p = pplay(i,k) call Psat_generic(v_t,v_p,metallicity,zpsat(i,k),zqs(i,k)) call Lcpdqsat_generic(v_t,v_p,zpsat(i,k),zqs(i,k),zdqs(i,k),zdlnpsat(i,k)) ENDDO ENDDO ! Calculate Gamma * Cp * dz: (gamma is the critical gradient) DO k = 2, nlayer DO i = 1, ngrid zdp = pplev(i,k)-pplev(i,k+1) zdpm = pplev(i,k-1)-pplev(i,k) ! gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp) & ! + RLVTT /(zdpm+zdp) * (zqs(i,k-1)*zdpm + zqs(i,k)*zdp) ) & !* (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) ) & ! / (1.0+ (zdqs(i,k-1)*zdpm + zdqs(i,k)*zdp)/(zdpm+zdp) ) ! general case where water is not a trace gas (JL13) v_zqs = (zqs(i,k-1)*zdpm + zqs(i,k)*zdp)/(zdpm+zdp) v_cptt2 = (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)/(zdpm+zdp) v_pratio = ((1.-zpsat(i,k-1)/pplay(i,k-1))*zdpm + (1.-zpsat(i,k)/pplay(i,k))*zdp)/(zdpm+zdp) v_dlnpsat = (zdlnpsat(i,k-1)*zdpm + zdlnpsat(i,k)*zdp)/(zdpm+zdp) gamcpdz(i,k) = ( (R/RCPD*v_cptt2*(1.- v_zqs) + RLVTT_generic*v_zqs) * (pplay(i,k-1)-pplay(i,k))/pplev(i,k) ) & / (((1.- v_zqs) + v_zqs * RCPV_generic/RCPD)*v_pratio + v_zqs * v_dlnpsat) ! Note that gamcpdz is defined as positive, so -gamcpdz is the real moist-adiabatic gradient [dT/dz]_ad ENDDO ENDDO ! calculate moist convection inhibition criterion IF (moist_convection_inhibition .and. (epsi_generic .gt. 1)) THEN ! GCS molecular weight is heavier than dry gas: ! inhibition of moist convection if vapor amount exceeds q_cri (Eq. 17 of Leconte et al. 2017) write(*,*) 'inhibition of moist convection if vapor amount exceeds a critical mixing ratio (see Leconte et al. 2017)' DO k = 1, nlayer DO i = 1, ngrid q_cri(i,k) = ( 1 / (1 - 1/epsi_generic)) * r * mugaz/1000. / delta_vapH * zt(i,k) ENDDO ENDDO ELSE ! GCS molecular weight is lighter than dry gas DO k = 1, nlayer DO i = 1, ngrid q_cri(i,k) = 2. ! vapor amount will never exceed q_cri, q_cri call becomes transparent in the next lines ENDDO ENDDO ENDIF !------------------------------------ defining the bottom (k1) and the top (k2) of the moist-convective column DO 9999 i = 1, ngrid itest(i) = .FALSE. ! print*,'we in the loop' ! stop k1 = 0 k2 = 1 810 CONTINUE ! look for k1, the base of the column k2 = k2 + 1 IF (k2 .GT. nlayer) GOTO 9999 zflo = v_cptt(i,k2-1) - v_cptt(i,k2) - gamcpdz(i,k2) zsat=(local_q(i,k2-1)-zqs(i,k2-1))*(pplev(i,k2-1)-pplev(i,k2)) & +(local_q(i,k2)-zqs(i,k2))*(pplev(i,k2)-pplev(i,k2+1)) ! if: (gradient is not steeper than moist-adiabat) or (level is not saturated) or (moist convection is inhibited because criterion is satisfied) ! then: no moist convection (GO TO 810) IF ( zflo.LE.0.0 .OR. zsat.LE.0.0 .OR. local_q(i,k2-1).GT.q_cri(i,k2-1)) GOTO 810 k1 = k2 - 1 itest(i) = .TRUE. 820 CONTINUE !! look for k2, the top of the column IF (k2 .EQ. nlayer) GOTO 821 k2p = k2 + 1 zsat=zsat+(pplev(i,k2p)-pplev(i,k2p+1))*(local_q(i,k2p)-zqs(i,k2p)) zflo = v_cptt(i,k2p-1) - v_cptt(i,k2p) - gamcpdz(i,k2p) IF (zflo.LE.0.0 .OR. zsat.LE.0.0 .OR. local_q(i,k2p).GT.q_cri(i,k2p)) GOTO 821 k2 = k2p GOTO 820 821 CONTINUE !------------------------------------------------------ local adjustment of the moist-convective column between k1 and k2 830 CONTINUE ! actual adjustment Do nn=1,niter v_cptj(k1) = 0.0 zdp = pplev(i,k1)-pplev(i,k1+1) v_cptjk1 = ( (1.0+zdqs(i,k1))*(v_cptt(i,k1)+v_cptj(k1)) + RLVTT_generic*(local_q(i,k1)-zqs(i,k1)) ) * zdp v_ssig = zdp * (1.0+zdqs(i,k1)) k1p = k1 + 1 DO k = k1p, k2 zdp = pplev(i,k)-pplev(i,k+1) v_cptj(k) = v_cptj(k-1) + gamcpdz(i,k) v_cptjk1 = v_cptjk1 + zdp * ( (1.0+zdqs(i, k))*(v_cptt(i,k)+v_cptj(k)) + RLVTT_generic*(local_q(i,k)-zqs(i,k)) ) v_ssig = v_ssig + zdp *(1.0+zdqs(i,k)) ENDDO ! this right here is where the adjustment is done??? DO k = k1, k2 cp_new_t(k) = v_cptjk1/v_ssig - v_cptj(k) cp_delta_t(k) = cp_new_t(k) - v_cptt(i,k) v_cptt(i,k)=cp_new_t(k) local_q(i,k) = zqs(i,k) + zdqs(i,k)*cp_delta_t(k)/RLVTT_generic local_t(i,k) = cp_new_t(k) / RCPD v_t = MAX(local_t(i,k),15.) v_p = pplay(i,k) call Psat_generic(v_t,v_p,metallicity,zpsat(i,k),zqs(i,k)) call Lcpdqsat_generic(v_t,v_p,zpsat(i,k),zqs(i,k),zdqs(i,k),zdlnpsat(i,k)) ENDDO Enddo ! nn=1,niter !--------------------------------------------------- sounding downwards ! -- we refine the prognostic variables in ! -- the layer about to be adjusted ! DO k = k1, k2 ! v_cptt(i,k) = RCPD * local_t(i,k) ! v_t = local_t(i,k) ! v_p = pplay(i,k) ! ! call Psat_water(v_t,v_p,zpsat,zqs(i,k)) ! call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k)) ! ENDDO DO k = 2, nlayer zdpm = pplev(i,k-1) - pplev(i,k) zdp = pplev(i,k) - pplev(i,k+1) ! gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp) & ! + RLVTT /(zdpm+zdp) * (zqs(i,k-1)*zdpm + zqs(i,k)*zdp) ) & ! * (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) ) & ! / (1.0+ (zdqs(i,k-1)*zdpm + zdqs(i,k)*zdp)/(zdpm+zdp) ) ! general case where water is not a trace gas (JL13) v_zqs = (zqs(i,k-1)*zdpm + zqs(i,k)*zdp)/(zdpm+zdp) v_cptt2 = (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)/(zdpm+zdp) v_pratio = ((1.-zpsat(i,k-1)/pplay(i,k-1))*zdpm + (1.-zpsat(i,k)/pplay(i,k))*zdp)/(zdpm+zdp) v_dlnpsat = (zdlnpsat(i,k-1)*zdpm + zdlnpsat(i,k)*zdp)/(zdpm+zdp) gamcpdz(i,k) = ( (R/RCPD*v_cptt2*(1.- v_zqs) + RLVTT_generic*v_zqs) * (pplay(i,k-1)-pplay(i,k))/pplev(i,k) ) & / (((1.- v_zqs) + v_zqs * RCPV_generic/RCPD)*v_pratio + v_zqs * v_dlnpsat) ENDDO ! Test to see if we've reached the bottom IF (k1 .EQ. 1) GOTO 841 ! yes we have! the bottom of moist-convective column is the bottom of the model zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1) zsat=(local_q(i,k1-1)-zqs(i,k1-1))*(pplev(i,k1-1)-pplev(i,k1)) & + (local_q(i,k1)-zqs(i,k1))*(pplev(i,k1)-pplev(i,k1+1)) IF (zflo.LE.0.0 .OR. zsat.LE.0.0 .OR. local_q(i,k1-1).GT.q_cri(i,k1-1)) GOTO 841 ! yes we have! ! the bottom of the moist-convective column is no more convective: GOTO 841 840 CONTINUE k1 = k1 - 1 IF (k1 .EQ. 1) GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995) ! bottom of moist-convective column is bottom of the model but is still convective: GOTO 830 zsat = zsat + (local_q(i,k1-1)-zqs(i,k1-1)) & *(pplev(i,k1-1)-pplev(i,k1)) zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1) IF (zflo.GT.0.0 .AND. zsat.GT.0.0 .AND. local_q(i,k1-1).LT.q_cri(i,k1-1)) THEN ! the bottom of the moist-convective column is still convective ! we continue to search for the non-convective bottom GOTO 840 ELSE ! the bottom of the moist-convective column is no more convective now, ! but, since the bottom of moist-convective column has changed, we must do again the adjustment GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995) ENDIF 841 CONTINUE GOTO 810 ! look for other layers higher up ! because there could be other moist-convective columns higher up (but separated from the one we have just calculated before) 9999 CONTINUE ! loop over all the points !----------------------------------------------------------------------- ! Determine the cloud fraction (hypothesis: nebulosity occurs ! where GCS vapor is reduced by adjustment): DO k = 1, nlayer DO i = 1, ngrid IF (itest(i)) THEN delta_q(i,k) = local_q(i,k) - zq(i,k) IF (delta_q(i,k).LT.0.) rneb(i,k) = 1.0 ENDIF ENDDO ENDDO ! Distribute GCS condensates into cloudy liquid/solid condensates (hypothesis: ! liquid/solid condensates are distributed to areas where GCS vapor ! decreases and are distributed in proportion to this decrease): DO i = 1, ngrid IF (itest(i)) THEN zq1(i) = 0.0 zq2(i) = 0.0 ENDIF ENDDO DO k = 1, nlayer DO i = 1, ngrid IF (itest(i)) THEN zdp = pplev(i,k)-pplev(i,k+1) zq1(i) = zq1(i) - delta_q(i,k) * zdp zq2(i) = zq2(i) - MIN(0.0, delta_q(i,k)) * zdp ENDIF ENDDO ENDDO DO k = 1, nlayer DO i = 1, ngrid IF (itest(i)) THEN IF (zq2(i).NE.0.0) d_ql(i,k) = - MIN(0.0,delta_q(i,k))*zq1(i)/zq2(i) ENDIF ENDDO ENDDO DO k = 1, nlayer DO i = 1, ngrid local_q(i, k) = MAX(local_q(i, k), seuil_vap) ENDDO ENDDO DO k = 1, nlayer DO i = 1, ngrid d_t(i,k) = local_t(i,k) - zt(i,k) d_q(i,k) = local_q(i,k) - zq(i,k) ENDDO ENDDO ! now subroutine -----> GCM variables DO k = 1, nlayer DO i = 1, ngrid pdtmana(i,k) = d_t(i,k)/ptimestep pdqmana(i,k,i_vap_generic) = d_q(i,k)/ptimestep pdqmana(i,k,i_ice_generic) = d_ql(i,k)/ptimestep ENDDO ENDDO end subroutine moistadj_generic