source: trunk/LMDZ.GENERIC/libf/phystd/moistadj_generic.F90 @ 3523

Last change on this file since 3523 was 3279, checked in by jleconte, 8 months ago

flag moist_convection_inhibition (Noe C)

File size: 15.4 KB
Line 
1subroutine moistadj_generic(ngrid, nlayer, nq, pt, pq, pdq, pplev, pplay, pdtmana, pdqmana, ptimestep, rneb)
2
3   use generic_cloud_common_h
4   use generic_tracer_index_mod, only: generic_tracer_index
5   use tracer_h
6   use ioipsl_getin_p_mod, only: getin_p !-> to get the metallicity
7   use comcstfi_mod, only: r, cpp, mugaz
8   use callkeys_mod, only: moist_convection_inhibition
9
10   implicit none
11
12
13!=====================================================================
14!     
15!     Purpose
16!     -------
17!     Calculates moist convective adjustment by the method of Manabe.
18!     
19!     Authors
20!     -------
21!     Adapted from the moistadj.F90 routine
22!     for generic condensable species (GCS) tracers
23!     by Noe CLEMENT (2023)
24!     
25!=====================================================================
26
27      INTEGER,INTENT(IN) :: ngrid, nlayer, nq
28
29      REAL,INTENT(IN) :: pt(ngrid,nlayer) ! temperature (K)
30      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracer (kg/kg)
31      REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq)
32      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
33      REAL,INTENT(IN) :: pplay(ngrid,nlayer)   ! mid-layer pressure (Pa)
34      REAL,INTENT(IN) :: ptimestep ! physics timestep (s)
35      REAL,INTENT(OUT) :: pdqmana(ngrid,nlayer,nq)  ! tracer tendencies (kg/kg.s-1)
36      REAL,INTENT(OUT) :: pdtmana(ngrid,nlayer) ! temperature increment(K/s)
37      REAL,INTENT(OUT) :: rneb(ngrid,nlayer) ! cloud fraction
38
39      ! Options :
40      !real, save :: metallicity ! metallicity of planet
41      REAL,SAVE :: metallicity = 0.0
42!$OMP THREADPRIVATE(metallicity)
43
44      ! local variables
45      REAL zt(ngrid,nlayer)      ! temperature (K)
46      REAL zq(ngrid,nlayer)      ! humidite specifique (kg/kg)
47
48      REAL d_t(ngrid,nlayer)     ! temperature increment
49      REAL d_q(ngrid,nlayer)     ! incrementation pour vapeur d'eau
50      REAL d_ql(ngrid,nlayer)    ! incrementation pour l'eau liquide
51
52      ! REAL t_coup
53      ! PARAMETER (t_coup=234.0)
54      REAL seuil_vap
55      PARAMETER (seuil_vap=1.0E-10)
56
57      ! Local variables
58
59      logical call_ice_vap_generic ! to call only one time the ice/vap pair of a tracer
60      integer igcm_generic_vap, igcm_generic_ice ! index of the vap and ice of generic_tracer
61
62      INTEGER i, k, iq, nn
63      INTEGER, PARAMETER :: niter = 1
64      INTEGER k1, k1p, k2, k2p
65      LOGICAL itest(ngrid)
66      REAL delta_q(ngrid, nlayer)
67      DOUBLE PRECISION :: cp_new_t(nlayer), v_cptt(ngrid,nlayer)
68      REAL cp_delta_t(nlayer)
69      DOUBLE PRECISION :: v_cptj(nlayer), v_cptjk1, v_ssig
70      REAL v_p, v_t, v_zqs,v_cptt2,v_pratio,v_dlnpsat
71      REAL zqs(ngrid,nlayer), zdqs(ngrid,nlayer),zpsat(ngrid,nlayer),zdlnpsat(ngrid,nlayer)
72      REAL zq1(ngrid), zq2(ngrid)
73      DOUBLE PRECISION :: gamcpdz(ngrid,2:nlayer)
74      DOUBLE PRECISION :: zdp, zdpm
75
76      real q_cri(ngrid,nlayer) ! moist convection inhibition criterion
77
78      REAL zsat ! super-saturation
79      REAL zflo ! flotabilite
80
81      DOUBLE PRECISION :: local_q(ngrid,nlayer),local_t(ngrid,nlayer)
82
83      REAL zdelta, zcor, zcvm5
84
85      REAL dEtot, dqtot, masse ! conservation diagnostics
86      real dL1tot, dL2tot
87
88!     Indices of generic vapour and ice tracers
89      real,save :: RCPD=0.0
90      INTEGER,SAVE :: i_vap_generic=0  ! Generic Condensable Species vapour
91      INTEGER,SAVE :: i_ice_generic=0  ! Generic Condensable Species ice
92!$OMP THREADPRIVATE(i_vap_generic,i_ice_generic,RCPD)
93
94      LOGICAL,SAVE :: firstcall=.TRUE.
95!$OMP THREADPRIVATE(firstcall)
96
97      IF (firstcall) THEN
98
99         RCPD = cpp
100
101         write(*,*) "value for metallicity? "
102         metallicity=0.0 ! default value
103         call getin_p("metallicity",metallicity)
104         write(*,*) " metallicity = ",metallicity
105
106         do iq=1, nq
107            call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,call_ice_vap_generic)
108            if(call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
109               write(*,*) "moistadjustment : taking back the values you have set in 'table_tracers_condensable' for thermodynamics. If you have forgotten one, it will crash."
110               m = constants_mass(iq)
111               delta_vapH = constants_delta_vapH(iq)
112               epsi_generic = constants_epsi_generic(iq)
113               RLVTT_generic = constants_RLVTT_generic(iq)
114               RCPV_generic = constants_RCPV_generic(iq)
115               Tref = constants_Tref(iq)
116
117               write(*,*) noms(igcm_generic_vap),", q_cri at ", Tref, "K (in kg/kg): ", ( 1 / (1 - 1/epsi_generic)) * (r * mugaz/1000.) / delta_vapH * Tref
118
119            endif
120         enddo
121         
122         i_vap_generic=igcm_generic_vap
123         i_ice_generic=igcm_generic_ice
124       
125         write(*,*) "rain: i_ice_generic=",i_ice_generic
126         write(*,*) "      i_vap_generic=",i_vap_generic
127
128         firstcall = .FALSE.
129      ENDIF
130
131      ! GCM -----> subroutine variables
132      zq(1:ngrid,1:nlayer)    = pq(1:ngrid,1:nlayer,i_vap_generic)+ pdq(1:ngrid,1:nlayer,i_vap_generic)*ptimestep
133      zt(1:ngrid,1:nlayer)    = pt(1:ngrid,1:nlayer)
134      pdqmana(1:ngrid,1:nlayer,1:nq)=0.0
135
136      DO k = 1, nlayer
137         DO i = 1, ngrid
138            if(zq(i,k).lt.0.)then
139               zq(i,k)=0.0
140            endif
141         ENDDO
142      ENDDO
143     
144      local_q(1:ngrid,1:nlayer) = zq(1:ngrid,1:nlayer)
145      local_t(1:ngrid,1:nlayer) = zt(1:ngrid,1:nlayer)
146      rneb(1:ngrid,1:nlayer) = 0.0
147      d_ql(1:ngrid,1:nlayer) = 0.0
148      d_t(1:ngrid,1:nlayer)  = 0.0
149      d_q(1:ngrid,1:nlayer)  = 0.0
150
151      !     Calculate v_cptt
152      DO k = 1, nlayer
153         DO i = 1, ngrid
154            v_cptt(i,k) = RCPD * local_t(i,k)
155            v_t = MAX(local_t(i,k),15.)
156            v_p = pplay(i,k)
157
158            call Psat_generic(v_t,v_p,metallicity,zpsat(i,k),zqs(i,k))
159            call Lcpdqsat_generic(v_t,v_p,zpsat(i,k),zqs(i,k),zdqs(i,k),zdlnpsat(i,k))
160
161         ENDDO
162      ENDDO
163
164      ! Calculate Gamma * Cp * dz: (gamma is the critical gradient)
165      DO k = 2, nlayer
166         DO i = 1, ngrid
167            zdp = pplev(i,k)-pplev(i,k+1)
168            zdpm = pplev(i,k-1)-pplev(i,k)
169            ! gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)          &
170            !     +  RLVTT /(zdpm+zdp) * (zqs(i,k-1)*zdpm + zqs(i,k)*zdp) )                           &
171            !* (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) )                                          &
172            !     / (1.0+ (zdqs(i,k-1)*zdpm + zdqs(i,k)*zdp)/(zdpm+zdp) )                   
173            ! general case where water is not a trace gas (JL13)
174            v_zqs     = (zqs(i,k-1)*zdpm + zqs(i,k)*zdp)/(zdpm+zdp)
175            v_cptt2   = (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)/(zdpm+zdp)
176            v_pratio  = ((1.-zpsat(i,k-1)/pplay(i,k-1))*zdpm + (1.-zpsat(i,k)/pplay(i,k))*zdp)/(zdpm+zdp)
177            v_dlnpsat = (zdlnpsat(i,k-1)*zdpm + zdlnpsat(i,k)*zdp)/(zdpm+zdp)
178            gamcpdz(i,k) = ( (R/RCPD*v_cptt2*(1.- v_zqs) + RLVTT_generic*v_zqs) * (pplay(i,k-1)-pplay(i,k))/pplev(i,k) )  &
179               / (((1.- v_zqs) + v_zqs * RCPV_generic/RCPD)*v_pratio + v_zqs  * v_dlnpsat)               
180            ! Note that gamcpdz is defined as positive, so -gamcpdz is the real moist-adiabatic gradient [dT/dz]_ad
181         ENDDO
182      ENDDO
183     
184      ! calculate moist convection inhibition criterion
185      IF (moist_convection_inhibition .and. (epsi_generic .gt. 1)) THEN ! GCS molecular weight is heavier than dry gas:
186         ! inhibition of moist convection if vapor amount exceeds q_cri (Eq. 17 of Leconte et al. 2017)
187         write(*,*) 'inhibition of moist convection if vapor amount exceeds a critical mixing ratio (see Leconte et al. 2017)'
188         DO k = 1, nlayer
189            DO i = 1, ngrid
190               q_cri(i,k) = ( 1 / (1 - 1/epsi_generic)) * r * mugaz/1000. / delta_vapH * zt(i,k)
191            ENDDO
192         ENDDO
193      ELSE ! GCS molecular weight is lighter than dry gas
194         DO k = 1, nlayer
195            DO i = 1, ngrid
196               q_cri(i,k) = 2. ! vapor amount will never exceed q_cri, q_cri call becomes transparent in the next lines
197            ENDDO
198         ENDDO
199      ENDIF
200
201
202!------------------------------------ defining the bottom (k1) and the top (k2) of the moist-convective column
203      DO 9999 i = 1, ngrid
204
205      itest(i) = .FALSE.
206
207      ! print*,'we in the loop'
208      ! stop   
209
210      k1 = 0
211      k2 = 1
212
213   810 CONTINUE ! look for k1, the base of the column
214      k2 = k2 + 1
215      IF (k2 .GT. nlayer) GOTO 9999
216      zflo = v_cptt(i,k2-1) - v_cptt(i,k2) - gamcpdz(i,k2)
217      zsat=(local_q(i,k2-1)-zqs(i,k2-1))*(pplev(i,k2-1)-pplev(i,k2))   &
218         +(local_q(i,k2)-zqs(i,k2))*(pplev(i,k2)-pplev(i,k2+1))
219
220      ! if: (gradient is not steeper than moist-adiabat) or (level is not saturated) or (moist convection is inhibited because criterion is satisfied)
221      ! then: no moist convection (GO TO 810)
222      IF ( zflo.LE.0.0 .OR. zsat.LE.0.0 .OR. local_q(i,k2-1).GT.q_cri(i,k2-1)) GOTO 810
223      k1 = k2 - 1
224      itest(i) = .TRUE.
225
226   820 CONTINUE !! look for k2, the top of the column
227      IF (k2 .EQ. nlayer) GOTO 821
228      k2p = k2 + 1
229      zsat=zsat+(pplev(i,k2p)-pplev(i,k2p+1))*(local_q(i,k2p)-zqs(i,k2p))
230      zflo = v_cptt(i,k2p-1) - v_cptt(i,k2p) - gamcpdz(i,k2p)
231
232      IF (zflo.LE.0.0 .OR. zsat.LE.0.0 .OR. local_q(i,k2p).GT.q_cri(i,k2p)) GOTO 821
233      k2 = k2p
234      GOTO 820
235   821 CONTINUE
236
237!------------------------------------------------------ local adjustment of the moist-convective column between k1 and k2
238   830 CONTINUE ! actual adjustment
239   Do nn=1,niter
240      v_cptj(k1) = 0.0
241      zdp = pplev(i,k1)-pplev(i,k1+1)
242      v_cptjk1 = ( (1.0+zdqs(i,k1))*(v_cptt(i,k1)+v_cptj(k1)) + RLVTT_generic*(local_q(i,k1)-zqs(i,k1)) ) * zdp
243      v_ssig = zdp * (1.0+zdqs(i,k1))
244
245      k1p = k1 + 1
246      DO k = k1p, k2
247         zdp = pplev(i,k)-pplev(i,k+1)
248         v_cptj(k) = v_cptj(k-1) + gamcpdz(i,k)
249         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)) )       
250         v_ssig = v_ssig + zdp *(1.0+zdqs(i,k))
251      ENDDO
252
253
254      ! this right here is where the adjustment is done???
255      DO k = k1, k2
256         cp_new_t(k) = v_cptjk1/v_ssig - v_cptj(k)
257         cp_delta_t(k) = cp_new_t(k) - v_cptt(i,k)
258         v_cptt(i,k)=cp_new_t(k)
259         local_q(i,k) = zqs(i,k) + zdqs(i,k)*cp_delta_t(k)/RLVTT_generic
260         local_t(i,k) = cp_new_t(k) / RCPD
261
262         v_t = MAX(local_t(i,k),15.)
263         v_p = pplay(i,k)
264
265         call Psat_generic(v_t,v_p,metallicity,zpsat(i,k),zqs(i,k))
266         call Lcpdqsat_generic(v_t,v_p,zpsat(i,k),zqs(i,k),zdqs(i,k),zdlnpsat(i,k))
267
268      ENDDO
269   Enddo ! nn=1,niter
270
271
272!--------------------------------------------------- sounding downwards
273!              -- we refine the prognostic variables in
274!              -- the layer about to be adjusted
275
276!      DO k = k1, k2
277!         v_cptt(i,k) = RCPD * local_t(i,k)
278!         v_t = local_t(i,k)
279!         v_p = pplay(i,k)
280!       
281!         call Psat_water(v_t,v_p,zpsat,zqs(i,k))
282!         call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k))
283!      ENDDO
284
285      DO k = 2, nlayer
286         zdpm = pplev(i,k-1) - pplev(i,k)
287         zdp = pplev(i,k) - pplev(i,k+1)
288         ! gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)          &
289         ! +  RLVTT /(zdpm+zdp) * (zqs(i,k-1)*zdpm + zqs(i,k)*zdp) )                           &
290         ! * (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) )                                          &
291         ! / (1.0+ (zdqs(i,k-1)*zdpm + zdqs(i,k)*zdp)/(zdpm+zdp) )                   
292         ! general case where water is not a trace gas (JL13)
293         v_zqs     = (zqs(i,k-1)*zdpm + zqs(i,k)*zdp)/(zdpm+zdp)
294         v_cptt2   = (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)/(zdpm+zdp)
295         v_pratio  = ((1.-zpsat(i,k-1)/pplay(i,k-1))*zdpm + (1.-zpsat(i,k)/pplay(i,k))*zdp)/(zdpm+zdp)
296         v_dlnpsat = (zdlnpsat(i,k-1)*zdpm + zdlnpsat(i,k)*zdp)/(zdpm+zdp)
297         gamcpdz(i,k) = ( (R/RCPD*v_cptt2*(1.- v_zqs) + RLVTT_generic*v_zqs) * (pplay(i,k-1)-pplay(i,k))/pplev(i,k) )  &
298               / (((1.- v_zqs) + v_zqs * RCPV_generic/RCPD)*v_pratio + v_zqs  * v_dlnpsat)               
299      ENDDO
300
301!     Test to see if we've reached the bottom
302
303      IF (k1 .EQ. 1) GOTO 841 ! yes we have! the bottom of moist-convective column is the bottom of the model
304
305      zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
306      zsat=(local_q(i,k1-1)-zqs(i,k1-1))*(pplev(i,k1-1)-pplev(i,k1))   &
307         + (local_q(i,k1)-zqs(i,k1))*(pplev(i,k1)-pplev(i,k1+1))
308
309      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!
310      ! the bottom of the moist-convective column is no more convective: GOTO 841
311
312   840 CONTINUE
313      k1 = k1 - 1
314      IF (k1 .EQ. 1) GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
315      ! bottom of moist-convective column is bottom of the model but is still convective: GOTO 830
316      zsat = zsat + (local_q(i,k1-1)-zqs(i,k1-1))               &
317                  *(pplev(i,k1-1)-pplev(i,k1))
318      zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
319      IF (zflo.GT.0.0 .AND. zsat.GT.0.0 .AND. local_q(i,k1-1).LT.q_cri(i,k1-1)) THEN
320         ! the bottom of the moist-convective column is still convective
321         ! we continue to search for the non-convective bottom
322         GOTO 840
323      ELSE
324         ! the bottom of the moist-convective column is no more convective now,
325         ! but, since the bottom of moist-convective column has changed, we must do again the adjustment
326         GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
327      ENDIF
328   841 CONTINUE
329
330      GOTO 810 ! look for other layers higher up
331      ! because there could be other moist-convective columns higher up (but separated from the one we have just calculated before)
332
3339999 CONTINUE ! loop over all the points
334
335!-----------------------------------------------------------------------
336! Determine the cloud fraction (hypothesis: nebulosity occurs
337! where GCS vapor is reduced by adjustment):
338
339      DO k = 1, nlayer
340      DO i = 1, ngrid
341         IF (itest(i)) THEN
342         delta_q(i,k) = local_q(i,k) - zq(i,k)
343         IF (delta_q(i,k).LT.0.) rneb(i,k)  = 1.0
344         ENDIF
345      ENDDO
346      ENDDO
347
348! Distribute GCS condensates into cloudy liquid/solid condensates (hypothesis:
349! liquid/solid condensates are distributed to areas where GCS vapor
350! decreases and are distributed in proportion to this decrease):
351      DO i = 1, ngrid
352         IF (itest(i)) THEN
353         zq1(i) = 0.0
354         zq2(i) = 0.0
355         ENDIF
356      ENDDO
357      DO k = 1, nlayer
358      DO i = 1, ngrid
359         IF (itest(i)) THEN
360         zdp = pplev(i,k)-pplev(i,k+1)
361         zq1(i) = zq1(i) - delta_q(i,k) * zdp
362         zq2(i) = zq2(i) - MIN(0.0, delta_q(i,k)) * zdp
363         ENDIF
364      ENDDO
365      ENDDO
366      DO k = 1, nlayer
367      DO i = 1, ngrid
368         IF (itest(i)) THEN
369            IF (zq2(i).NE.0.0) d_ql(i,k) = - MIN(0.0,delta_q(i,k))*zq1(i)/zq2(i)
370         ENDIF
371      ENDDO
372      ENDDO
373
374      DO k = 1, nlayer
375      DO i = 1, ngrid
376         local_q(i, k) = MAX(local_q(i, k), seuil_vap)
377      ENDDO
378      ENDDO
379
380      DO k = 1, nlayer
381      DO i = 1, ngrid
382         d_t(i,k) = local_t(i,k) - zt(i,k)
383         d_q(i,k) = local_q(i,k) - zq(i,k)
384      ENDDO
385      ENDDO
386
387      ! now subroutine -----> GCM variables
388      DO k = 1, nlayer
389         DO i = 1, ngrid
390           
391            pdtmana(i,k)       = d_t(i,k)/ptimestep
392            pdqmana(i,k,i_vap_generic) = d_q(i,k)/ptimestep
393            pdqmana(i,k,i_ice_generic) = d_ql(i,k)/ptimestep
394         
395         ENDDO
396      ENDDO
397
398end subroutine moistadj_generic
Note: See TracBrowser for help on using the repository browser.