source: LMDZ6/branches/cirrus/libf/phylmdiso/climb_hq_mod.F90 @ 5405

Last change on this file since 5405 was 4143, checked in by dcugnet, 3 years ago
  • Some variables are renamed or replaced by direct equivalents:
    • iso_indnum -> tracers(:)%iso_iName
    • niso_possibles -> niso
    • iqiso -> iqIsoPha ; index_trac -> itZonIso
    • ok_iso_verif -> isoCheck
    • ntraceurs_zone -> nzone ; ntraciso -> ntiso
    • qperemin -> min_qparent ; masseqmin -> min_qmass ; ratiomin -> min_ratio
  • Some renamed variables are only aliased with the older name (using USE <module>, ONLY: <oldName> => <newName>) in routines where they are repeated many times.
  • Few hard-coded indexes are now computed (examples: ilic, iso, ivap, irneb, iq_vap, iq_liq, iso_H2O, iso_HDO, iso_HTO, iso_O17, iso_O18).
  • The IF(isoCheck) test is now embedded in the check_isotopes_seq and check_isotopes_loc routines (lighter calling).
  • Property svn:keywords set to Id
File size: 28.1 KB
Line 
1MODULE climb_hq_mod
2!
3! Module to solve the verctical diffusion of "q" and "H";
4! specific humidity and potential energi.
5!
6  USE dimphy
7#ifdef ISO
8  USE infotrac_phy, ONLY: ntraciso=>ntiso ! ajout C Risi pour isos     
9#endif
10
11  IMPLICIT NONE
12  SAVE
13  PRIVATE
14  PUBLIC :: climb_hq_down, climb_hq_up, d_h_col_vdf, f_h_bnd
15
16  REAL, DIMENSION(:,:), ALLOCATABLE :: gamaq, gamah
17  !$OMP THREADPRIVATE(gamaq,gamah)
18  REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_Q, Dcoef_Q
19  !$OMP THREADPRIVATE(Ccoef_Q, Dcoef_Q)
20  REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_H, Dcoef_H
21  !$OMP THREADPRIVATE(Ccoef_H, Dcoef_H)
22  REAL, DIMENSION(:), ALLOCATABLE   :: Acoef_Q, Bcoef_Q
23  !$OMP THREADPRIVATE(Acoef_Q, Bcoef_Q)
24  REAL, DIMENSION(:), ALLOCATABLE   :: Acoef_H, Bcoef_H
25  !$OMP THREADPRIVATE(Acoef_H, Bcoef_H)
26  REAL, DIMENSION(:,:), ALLOCATABLE :: Kcoefhq
27  !$OMP THREADPRIVATE(Kcoefhq)
28  REAL, SAVE, DIMENSION(:,:), ALLOCATABLE :: h_old ! for diagnostics, h before solving diffusion
29  !$OMP THREADPRIVATE(h_old)
30  REAL, SAVE, DIMENSION(:), ALLOCATABLE :: d_h_col_vdf ! for diagnostics, vertical integral of enthalpy change
31  !$OMP THREADPRIVATE(d_h_col_vdf)
32  REAL, SAVE, DIMENSION(:), ALLOCATABLE :: f_h_bnd ! for diagnostics, enthalpy flux at surface
33  !$OMP THREADPRIVATE(f_h_bnd)
34#ifdef ISO
35  REAL, DIMENSION(:,:,:), ALLOCATABLE :: gamaxt
36  !$OMP THREADPRIVATE(gamaxt)
37  REAL, DIMENSION(:,:,:), ALLOCATABLE :: Ccoef_XT, Dcoef_XT
38  !$OMP THREADPRIVATE(Ccoef_XT, Dcoef_XT)
39  REAL, DIMENSION(:,:), ALLOCATABLE   :: Acoef_XT, Bcoef_XT
40  !$OMP THREADPRIVATE(Acoef_XT, Bcoef_XT)
41#endif
42
43CONTAINS
44!
45!****************************************************************************************
46!
47  SUBROUTINE climb_hq_down(knon, coefhq, paprs, pplay, &
48       delp, temp, q, dtime, &
49!!! nrlmd le 02/05/2011
50       Ccoef_H_out, Ccoef_Q_out, Dcoef_H_out, Dcoef_Q_out, &
51       Kcoef_hq_out, gama_q_out, gama_h_out, &
52!!!
53       Acoef_H_out, Acoef_Q_out, Bcoef_H_out, Bcoef_Q_out &
54#ifdef ISO
55            ,xt,  &
56       Ccoef_XT_out,Dcoef_XT_out,gama_xt_out,  &
57       Acoef_XT_out, Bcoef_XT_out &
58#endif               
59            )
60#ifdef ISOVERIF
61USE isotopes_mod, ONLY: iso_eau,iso_HDO
62!USE isotopes_verif_mod, ONLY: errmax, errmaxrel
63USE isotopes_verif_mod
64#endif
65
66! This routine calculates recursivly the coefficients C and D
67! for the quantity X=[Q,H] in equation X(k) = C(k) + D(k)*X(k-1), where k is
68! the index of the vertical layer.
69!
70! Input arguments
71!****************************************************************************************
72    INTEGER, INTENT(IN)                      :: knon
73    REAL, DIMENSION(klon,klev), INTENT(IN)   :: coefhq
74    REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay
75    REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
76    REAL, DIMENSION(klon,klev), INTENT(IN)   :: temp, delp  ! temperature
77    REAL, DIMENSION(klon,klev), INTENT(IN)   :: q
78    REAL, INTENT(IN)                         :: dtime
79#ifdef ISO
80    REAL, DIMENSION(ntraciso,klon,klev), INTENT(IN)   :: xt
81#endif
82
83
84! Output arguments
85!****************************************************************************************
86    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef_H_out
87    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef_Q_out
88    REAL, DIMENSION(klon), INTENT(OUT)       :: Bcoef_H_out
89    REAL, DIMENSION(klon), INTENT(OUT)       :: Bcoef_Q_out
90#ifdef ISO
91    REAL, DIMENSION(ntraciso,klon), INTENT(OUT)       :: Acoef_XT_out
92    REAL, DIMENSION(ntraciso,klon), INTENT(OUT)       :: Bcoef_XT_out
93#endif
94
95!!! nrlmd le 02/05/2011
96    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Ccoef_H_out
97    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Ccoef_Q_out
98    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Dcoef_H_out
99    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Dcoef_Q_out
100    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Kcoef_hq_out
101    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: gama_q_out
102    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: gama_h_out
103#ifdef ISO
104    REAL, DIMENSION(ntraciso,klon,klev), INTENT(OUT)  :: Ccoef_XT_out
105    REAL, DIMENSION(ntraciso,klon,klev), INTENT(OUT)  :: Dcoef_XT_out
106    REAL, DIMENSION(ntraciso,klon,klev), INTENT(OUT)  :: gama_xt_out
107#endif
108!!!
109
110! Local variables
111!****************************************************************************************
112    LOGICAL, SAVE                            :: first=.TRUE.
113    !$OMP THREADPRIVATE(first)
114! JLD now renamed h_old and declared in module
115!    REAL, DIMENSION(klon,klev)               :: local_H
116    REAL, DIMENSION(klon)                    :: psref
117    REAL                                     :: delz, pkh
118    INTEGER                                  :: k, i, ierr
119
120#ifdef ISO
121    real, DIMENSION(klon,2:klev) :: gamaxt_tmp
122    real, DIMENSION(klon,klev) :: xt_tmp
123    real, DIMENSION(klon,klev) ::  Ccoef_XT_tmp,Dcoef_XT_tmp
124    real, DIMENSION(klon) ::  Acoef_XT_tmp,Bcoef_XT_tmp
125    integer ixt
126#endif
127! Include
128!****************************************************************************************
129    INCLUDE "YOMCST.h"
130    INCLUDE "compbl.h"   
131   
132#ifdef ISO
133#ifdef ISOVERIF
134   if (iso_eau.gt.0) then
135       do k = 1, klev
136          DO i = 1, knon
137            call iso_verif_egalite_choix( &
138                xt(iso_eau,i,k),q(i,k), &
139                'climb_hq 100',errmax,errmaxrel)
140          enddo
141       enddo
142   endif ! if (iso_eau.gt.0) then
143#endif
144#endif
145
146!****************************************************************************************
147! 1)
148! Allocation at first time step only
149!   
150!****************************************************************************************
151
152    IF (first) THEN
153       first=.FALSE.
154       ALLOCATE(Ccoef_Q(klon,klev), STAT=ierr)
155       IF ( ierr /= 0 )  PRINT*,' pb in allloc Ccoef_Q, ierr=', ierr
156       
157       ALLOCATE(Dcoef_Q(klon,klev), STAT=ierr)
158       IF ( ierr /= 0 )  PRINT*,' pb in allloc Dcoef_Q, ierr=', ierr
159       
160       ALLOCATE(Ccoef_H(klon,klev), STAT=ierr)
161       IF ( ierr /= 0 )  PRINT*,' pb in allloc Ccoef_H, ierr=', ierr
162       
163       ALLOCATE(Dcoef_H(klon,klev), STAT=ierr)
164       IF ( ierr /= 0 )  PRINT*,' pb in allloc Dcoef_H, ierr=', ierr
165       
166       ALLOCATE(Acoef_Q(klon), Bcoef_Q(klon), Acoef_H(klon), Bcoef_H(klon), STAT=ierr)
167       IF ( ierr /= 0 )  PRINT*,' pb in allloc Acoef_X and Bcoef_X, ierr=', ierr
168       
169       ALLOCATE(Kcoefhq(klon,klev), STAT=ierr)
170       IF ( ierr /= 0 )  PRINT*,' pb in allloc Kcoefhq, ierr=', ierr
171       
172       ALLOCATE(gamaq(1:klon,2:klev), STAT=ierr)
173       IF ( ierr /= 0 ) PRINT*,' pb in allloc gamaq, ierr=', ierr
174       
175       ALLOCATE(gamah(1:klon,2:klev), STAT=ierr)
176       IF ( ierr /= 0 ) PRINT*,' pb in allloc gamah, ierr=', ierr
177
178#ifdef ISO
179       ALLOCATE(Ccoef_XT(ntraciso,klon,klev), STAT=ierr)
180       IF ( ierr /= 0 )  PRINT*,' pb in allloc Ccoef_XT, ierr=', ierr
181       
182       ALLOCATE(Dcoef_XT(ntraciso,klon,klev), STAT=ierr)
183       IF ( ierr /= 0 )  PRINT*,' pb in allloc Dcoef_XT, ierr=', ierr
184       
185       ALLOCATE(Acoef_XT(ntraciso,klon), Bcoef_XT(ntraciso,klon), STAT=ierr)
186       IF ( ierr /= 0 )  PRINT*,' pb in allloc Acoef_XT and Bcoef_XT, ierr=', ierr
187       
188        ALLOCATE(gamaxt(ntraciso,1:klon,2:klev), STAT=ierr)
189       IF ( ierr /= 0 ) PRINT*,' pb in allloc gamaxt, ierr=', ierr
190#endif
191       
192       ALLOCATE(h_old(klon,klev), STAT=ierr)
193       IF ( ierr /= 0 )  PRINT*,' pb in allloc h_old, ierr=', ierr
194       
195       ALLOCATE(d_h_col_vdf(klon), STAT=ierr)
196       IF ( ierr /= 0 )  PRINT*,' pb in allloc d_h_col_vdf, ierr=', ierr
197       
198       ALLOCATE(f_h_bnd(klon), STAT=ierr)
199       IF ( ierr /= 0 )  PRINT*,' pb in allloc f_h_bnd, ierr=', ierr
200    END IF
201
202!****************************************************************************************
203! 2)
204! Definition of the coeficient K
205!
206!****************************************************************************************
207    Kcoefhq(:,:) = 0.0
208    DO k = 2, klev
209       DO i = 1, knon
210          Kcoefhq(i,k) = &
211               coefhq(i,k)*RG*RG*dtime /(pplay(i,k-1)-pplay(i,k)) &
212               *(paprs(i,k)*2/(temp(i,k)+temp(i,k-1))/RD)**2
213       ENDDO
214    ENDDO
215
216!****************************************************************************************
217! 3)
218! Calculation of gama for "Q" and "H"
219!
220!****************************************************************************************
221!   surface pressure is used as reference
222    psref(:) = paprs(:,1)
223
224!   definition of gama
225    IF (iflag_pbl == 1) THEN
226       gamaq(:,:) = 0.0
227       gamah(:,:) = -1.0e-03
228       gamah(:,2) = -2.5e-03
229#ifdef ISO
230       do ixt=1,ntraciso
231        gamaxt(:,:,:) = 0.0
232       enddo ! do ixt=1,ntraciso
233#endif
234 
235! conversion de gama
236       DO k = 2, klev
237          DO i = 1, knon
238             delz = RD * (temp(i,k-1)+temp(i,k)) / &
239                    2.0 / RG / paprs(i,k) * (pplay(i,k-1)-pplay(i,k))
240             pkh  = (psref(i)/paprs(i,k))**RKAPPA
241         
242! convertie gradient verticale d'humidite specifique en difference d'humidite specifique entre centre de couches
243             gamaq(i,k) = gamaq(i,k) * delz   
244! convertie gradient verticale de temperature en difference de temperature potentielle entre centre de couches
245             gamah(i,k) = gamah(i,k) * delz * RCPD * pkh
246#ifdef ISO
247             do ixt=1,ntraciso
248              gamaxt(ixt,i,k) = gamaxt(ixt,i,k) * delz
249             enddo
250#endif
251          ENDDO
252       ENDDO
253
254    ELSE
255       gamaq(:,:) = 0.0
256       gamah(:,:) = 0.0
257#ifdef ISO
258       do ixt=1,ntraciso
259        gamaxt(:,:,:) = 0.0
260       enddo ! do ixt=1,ntraciso
261#endif
262    ENDIF
263   
264#ifdef ISO
265#ifdef ISOVERIF
266        do k = 2, klev
267          DO i = 1, knon
268            call iso_verif_egalite_choix( &
269                gamaxt(iso_eau,i,k),gamaq(i,k), &
270                'climb_hq 209',errmax,errmaxrel)
271          enddo
272       enddo
273#endif
274#endif
275   
276
277!****************************************************************************************   
278! 4)
279! Calculte the coefficients C and D for specific humidity, q
280!
281!****************************************************************************************
282   
283    CALL calc_coef(knon, Kcoefhq(:,:), gamaq(:,:), delp(:,:), q(:,:), &
284         Ccoef_Q(:,:), Dcoef_Q(:,:), Acoef_Q, Bcoef_Q)
285
286
287#ifdef ISO
288        do ixt=1,ntraciso       
289        ! compression
290        do k = 2, klev
291          DO i = 1, knon
292            gamaxt_tmp(i,k)=gamaxt(ixt,i,k)
293          enddo
294        enddo !do k = 2, klev
295        do k = 1, klev
296          DO i = 1, knon
297            xt_tmp(i,k)=xt(ixt,i,k)
298          enddo
299        enddo !do k = 2, klev
300        !appel routine generique
301          CALL calc_coef(knon, Kcoefhq(:,:), gamaxt_tmp(:,:), delp(:,:), xt_tmp(:,:), &
302                Ccoef_XT_tmp(:,:), Dcoef_XT_tmp(:,:), Acoef_XT_tmp, Bcoef_XT_tmp) 
303
304        ! decompression
305        do k = 1, klev
306          DO i = 1, knon
307            Ccoef_XT(ixt,i,k)=Ccoef_XT_tmp(i,k)
308            Dcoef_XT(ixt,i,k)=Dcoef_XT_tmp(i,k)
309          enddo
310        enddo !do k = 2, klev
311        DO i = 1, knon
312            Acoef_XT(ixt,i)=Acoef_XT_tmp(i)
313            Bcoef_XT(ixt,i)=Bcoef_XT_tmp(i)
314        enddo
315       enddo ! do ixt=1,ntraciso
316#ifdef ISOVERIF
317        if (iso_eau.gt.0) then
318         do k = 1, klev
319          DO i = 1, knon
320            call iso_verif_egalite_choix( &
321                Ccoef_XT(iso_eau,i,k),Ccoef_Q(i,k), &
322                        'climb_hq 234c',errmax,errmaxrel)
323            call iso_verif_egalite_choix( &
324                Dcoef_XT(iso_eau,i,k),Dcoef_Q(i,k), &
325                        'climb_hq 234d',errmax,errmaxrel)
326          enddo !DO i = 1, knon
327         enddo !do k = 2, klev
328         DO i = 1, knon
329            call iso_verif_egalite_choix( &
330                Acoef_XT(iso_eau,i),Acoef_Q(i), &
331                        'climb_hq 234a',errmax,errmaxrel)
332            call iso_verif_egalite_choix( &
333                Bcoef_XT(iso_eau,i),Bcoef_Q(i), &
334                        'climb_hq 234b',errmax,errmaxrel)
335         enddo !DO i = 1, knon
336        endif !if (iso_eau.gt.0) then
337#endif
338#endif 
339!****************************************************************************************
340! 5)
341! Calculte the coefficients C and D for potentiel entalpie, H
342!
343!****************************************************************************************
344    h_old(:,:) = 0.0
345
346    DO k=1,klev
347       DO i = 1, knon
348          ! convertie la temperature en entalpie potentielle
349          h_old(i,k) = RCPD * temp(i,k) * &
350               (psref(i)/pplay(i,k))**RKAPPA
351       ENDDO
352    ENDDO
353
354    CALL calc_coef(knon, Kcoefhq(:,:), gamah(:,:), delp(:,:), h_old(:,:), &
355         Ccoef_H(:,:), Dcoef_H(:,:), Acoef_H, Bcoef_H)
356 
357!****************************************************************************************
358! 6)
359! Return the first layer in output variables
360!
361!****************************************************************************************
362    Acoef_H_out = Acoef_H
363    Bcoef_H_out = Bcoef_H
364    Acoef_Q_out = Acoef_Q
365    Bcoef_Q_out = Bcoef_Q
366#ifdef ISO
367    Acoef_XT_out = Acoef_XT
368    Bcoef_XT_out = Bcoef_XT
369#endif
370
371!****************************************************************************************
372! 7)
373! If Pbl is split, return also the other layers in output variables
374!
375!****************************************************************************************
376!!! jyg le 07/02/2012
377!!jyg       IF (mod(iflag_pbl_split,2) .eq.1) THEN
378       IF (mod(iflag_pbl_split,10) .ge.1) THEN
379!!! nrlmd le 02/05/2011
380    DO k= 1, klev
381      DO i= 1, klon
382        Ccoef_H_out(i,k) = Ccoef_H(i,k)
383        Dcoef_H_out(i,k) = Dcoef_H(i,k)
384        Ccoef_Q_out(i,k) = Ccoef_Q(i,k)
385        Dcoef_Q_out(i,k) = Dcoef_Q(i,k)
386        Kcoef_hq_out(i,k) = Kcoefhq(i,k)
387#ifdef ISO
388        do ixt=1,ntraciso
389          Ccoef_XT_out(ixt,i,k) = Ccoef_XT(ixt,i,k)
390          Dcoef_XT_out(ixt,i,k) = Dcoef_XT(ixt,i,k)             
391        enddo   
392#endif
393          IF (k.eq.1) THEN
394            gama_h_out(i,k)  = 0.
395            gama_q_out(i,k)  = 0.
396#ifdef ISO
397            do ixt=1,ntraciso
398              gama_xt_out(ixt,i,k)  = 0.
399            enddo   
400#endif
401          ELSE
402            gama_h_out(i,k)  = gamah(i,k)
403            gama_q_out(i,k)  = gamaq(i,k)
404#ifdef ISO
405            do ixt=1,ntraciso
406              gama_xt_out(ixt,i,k)  = gamaxt(ixt,i,k)
407            enddo   
408#endif
409          ENDIF
410      ENDDO
411    ENDDO
412!!!     
413       ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
414!!!
415
416  END SUBROUTINE climb_hq_down
417!
418!****************************************************************************************
419!
420  SUBROUTINE calc_coef(knon, Kcoef, gama, delp, X, Ccoef, Dcoef, Acoef, Bcoef)
421!
422! Calculate the coefficients C and D in : X(k) = C(k) + D(k)*X(k-1)
423! where X is H or Q, and k the vertical level k=1,klev
424!
425    INCLUDE "YOMCST.h"
426! Input arguments
427!****************************************************************************************
428    INTEGER, INTENT(IN)                      :: knon
429    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Kcoef, delp
430    REAL, DIMENSION(klon,klev), INTENT(IN)   :: X
431    REAL, DIMENSION(klon,2:klev), INTENT(IN) :: gama
432
433! Output arguments
434!****************************************************************************************
435    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef, Bcoef
436    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Ccoef, Dcoef
437
438! Local variables
439!****************************************************************************************
440    INTEGER                                  :: k, i
441    REAL                                     :: buf
442
443!****************************************************************************************
444! Niveau au sommet, k=klev
445!
446!****************************************************************************************
447    Ccoef(:,:) = 0.0
448    Dcoef(:,:) = 0.0
449
450    DO i = 1, knon
451       buf = delp(i,klev) + Kcoef(i,klev)
452       
453       Ccoef(i,klev) = (X(i,klev)*delp(i,klev) - Kcoef(i,klev)*gama(i,klev))/buf
454       Dcoef(i,klev) = Kcoef(i,klev)/buf
455    END DO
456
457
458!****************************************************************************************
459! Niveau  (klev-1) <= k <= 2
460!
461!****************************************************************************************
462
463    DO k=(klev-1),2,-1
464       DO i = 1, knon
465          buf = delp(i,k) + Kcoef(i,k) + Kcoef(i,k+1)*(1.-Dcoef(i,k+1))
466          Ccoef(i,k) = (X(i,k)*delp(i,k) + Kcoef(i,k+1)*Ccoef(i,k+1) + &
467               Kcoef(i,k+1)*gama(i,k+1) - Kcoef(i,k)*gama(i,k))/buf
468          Dcoef(i,k) = Kcoef(i,k)/buf
469       END DO
470    END DO
471
472!****************************************************************************************
473! Niveau k=1
474!
475!****************************************************************************************
476
477    DO i = 1, knon
478       buf = delp(i,1) + Kcoef(i,2)*(1.-Dcoef(i,2))
479       Acoef(i) = (X(i,1)*delp(i,1) + Kcoef(i,2)*(gama(i,2)+Ccoef(i,2)))/buf
480       Bcoef(i) = -1. * RG / buf
481    END DO
482
483  END SUBROUTINE calc_coef
484!
485!****************************************************************************************
486!
487  SUBROUTINE climb_hq_up(knon, dtime, t_old, q_old, &
488       flx_q1, flx_h1, paprs, pplay, &
489!!! nrlmd le 02/05/2011
490       Acoef_H_in, Acoef_Q_in, Bcoef_H_in, Bcoef_Q_in, &
491       Ccoef_H_in, Ccoef_Q_in, Dcoef_H_in, Dcoef_Q_in, &
492       Kcoef_hq_in, gama_q_in, gama_h_in, &
493!!!
494       flux_q, flux_h, d_q, d_t &
495#ifdef ISO
496       ,xt_old, flx_xt1, &
497       Acoef_XT_in,Bcoef_XT_in,Ccoef_XT_in,Dcoef_XT_in,gama_xt_in,  &
498       flux_xt, d_xt &
499#endif       
500       )
501
502#ifdef ISOVERIF
503USE infotrac_phy, ONLY: nzone
504USE isotopes_mod, ONLY: iso_eau,iso_HDO,iso_O18, ridicule
505USE isotopes_verif_mod
506#endif
507!
508! This routine calculates the flux and tendency of the specific humidity q and
509! the potential engergi H.
510! The quantities q and H are calculated according to
511! X(k) = C(k) + D(k)*X(k-1) for X=[q,H], where the coefficients
512! C and D are known from before and k is index of the vertical layer.
513!   
514
515! Input arguments
516!****************************************************************************************
517    INTEGER, INTENT(IN)                      :: knon
518    REAL, INTENT(IN)                         :: dtime
519    REAL, DIMENSION(klon,klev), INTENT(IN)   :: t_old, q_old
520    REAL, DIMENSION(klon), INTENT(IN)        :: flx_q1, flx_h1
521    REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
522    REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay
523
524!!! nrlmd le 02/05/2011
525    REAL, DIMENSION(klon), INTENT(IN)        :: Acoef_H_in,Acoef_Q_in, Bcoef_H_in, Bcoef_Q_in
526    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Ccoef_H_in, Ccoef_Q_in, Dcoef_H_in, Dcoef_Q_in
527    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Kcoef_hq_in, gama_q_in, gama_h_in
528#ifdef ISO
529    REAL, DIMENSION(ntraciso,klon), INTENT(IN)        :: Acoef_XT_in,  Bcoef_XT_in
530    REAL, DIMENSION(ntraciso,klon,klev), INTENT(IN)   :: Ccoef_XT_in,  Dcoef_XT_in
531    REAL, DIMENSION(ntraciso,klon,klev), INTENT(IN)   :: gama_xt_in
532#endif
533!!!
534
535! Output arguments
536!****************************************************************************************
537    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: flux_q, flux_h, d_q, d_t
538
539! Local variables
540!****************************************************************************************
541    LOGICAL, SAVE                            :: last=.FALSE.
542    REAL, DIMENSION(klon,klev)               :: h_new, q_new
543    REAL, DIMENSION(klon)                    :: psref         
544    INTEGER                                  :: k, i, ierr
545 
546#ifdef ISO
547    REAL, DIMENSION(ntraciso,klon,klev), INTENT(IN)   ::  xt_old
548    REAL, DIMENSION(ntraciso,klon,klev), INTENT(OUT)  :: flux_xt, d_xt
549    REAL, DIMENSION(ntraciso,klon,klev)               :: xt_new
550    REAL, DIMENSION(ntraciso,klon), INTENT(IN)        :: flx_xt1
551    integer ixt
552!#ifdef ISOVERIF
553!    integer iso_verif_noNaN_nostop
554!#endif
555#endif
556 
557! Include
558!****************************************************************************************
559    INCLUDE "YOMCST.h"
560    INCLUDE "compbl.h"       
561
562!****************************************************************************************
563! 1)
564! Definition of some variables
565    REAL, DIMENSION(klon,klev)               :: d_h, zairm
566!
567!****************************************************************************************
568
569
570#ifdef ISO
571#ifdef ISOVERIF
572       DO k = 1, klev
573        DO i = 1, knon 
574         if (iso_eau.gt.0) then
575           call iso_verif_egalite(xt_old(iso_eau,i,k), &
576    &          q_old(i,k),'climb_hq_mod 421')
577         endif
578#ifdef ISOTRAC
579         IF(nzone > 0) CALL iso_verif_traceur(xt_old(1,i,k),'climb_hq_mod 422')
580#endif
581        enddo
582       enddo
583#endif 
584#endif
585
586    flux_q(:,:) = 0.0
587    flux_h(:,:) = 0.0
588    d_q(:,:)    = 0.0
589    d_t(:,:)    = 0.0
590    d_h(:,:)    = 0.0
591    f_h_bnd(:)= 0.0
592#ifdef ISO
593    q_new(:,:)    = 0.0
594    flux_xt(:,:,:) = 0.0
595    d_xt(:,:,:)    = 0.0   
596    xt_new(:,:,:)    = 0.0
597#endif
598    psref(1:knon) = paprs(1:knon,1) 
599
600!!! jyg le 07/02/2012
601!!jyg       IF (mod(iflag_pbl_split,2) .eq.1) THEN
602       IF (mod(iflag_pbl_split,10) .ge.1) THEN
603!!! nrlmd le 02/05/2011
604    DO i = 1, knon
605      Acoef_H(i)=Acoef_H_in(i)
606      Acoef_Q(i)=Acoef_Q_in(i)
607      Bcoef_H(i)=Bcoef_H_in(i)
608      Bcoef_Q(i)=Bcoef_Q_in(i)
609#ifdef ISO
610      do ixt=1,ntraciso
611        Acoef_XT(ixt,i)=Acoef_XT_in(ixt,i)
612        Bcoef_XT(ixt,i)=Bcoef_XT_in(ixt,i)
613      enddo   
614#endif
615    ENDDO
616    DO k = 1, klev
617      DO i = 1, knon
618        Ccoef_H(i,k)=Ccoef_H_in(i,k)
619        Ccoef_Q(i,k)=Ccoef_Q_in(i,k)
620        Dcoef_H(i,k)=Dcoef_H_in(i,k)
621        Dcoef_Q(i,k)=Dcoef_Q_in(i,k)
622        Kcoefhq(i,k)=Kcoef_hq_in(i,k)
623#ifdef ISO
624      do ixt=1,ntraciso
625        Ccoef_XT(ixt,i,k)=Ccoef_XT_in(ixt,i,k)
626        Dcoef_XT(ixt,i,k)=Dcoef_XT_in(ixt,i,k)
627      enddo   
628#endif
629          IF (k.gt.1) THEN
630            gamah(i,k)=gama_h_in(i,k)
631            gamaq(i,k)=gama_q_in(i,k)
632#ifdef ISO
633            do ixt=1,ntraciso
634              gamaxt(ixt,i,k)=gama_xt_in(ixt,i,k)
635            enddo
636#endif
637          ENDIF
638      ENDDO
639    ENDDO
640!!!     
641       ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
642!!!
643
644!****************************************************************************************
645! 2)
646! Calculation of Q and H
647!
648!****************************************************************************************
649
650!- First layer
651    q_new(1:knon,1) = Acoef_Q(1:knon) + Bcoef_Q(1:knon)*flx_q1(1:knon)*dtime
652    h_new(1:knon,1) = Acoef_H(1:knon) + Bcoef_H(1:knon)*flx_h1(1:knon)*dtime
653    f_h_bnd(1:knon) = flx_h1(1:knon)
654#ifdef ISO
655    do ixt=1,ntraciso
656      xt_new(ixt,1:knon,1) = Acoef_XT(ixt,1:knon) + Bcoef_XT(ixt,1:knon)*flx_xt1(ixt,1:knon)*dtime
657    enddo ! do ixt=1,ntraciso
658#endif
659!- All the other layers
660    DO k = 2, klev
661       DO i = 1, knon
662          q_new(i,k) = Ccoef_Q(i,k) + Dcoef_Q(i,k)*q_new(i,k-1)
663          h_new(i,k) = Ccoef_H(i,k) + Dcoef_H(i,k)*h_new(i,k-1)
664#ifdef ISO
665        do ixt=1,ntraciso
666          xt_new(ixt,i,k) = Ccoef_XT(ixt,i,k) + Dcoef_XT(ixt,i,k)*xt_new(ixt,i,k-1)
667        enddo ! do ixt=1,ntraciso
668#endif
669       END DO
670    END DO
671
672#ifdef ISO
673#ifdef ISOVERIF
674    DO k = 1, klev
675       DO i = 1, knon     
676        do ixt=1,ntraciso
677         if (iso_verif_noNaN_nostop(xt_new(ixt,i,k),'climb_hq 507').eq.1) then
678             write(*,*) 'Acoef_XT(ixt,i)=',Acoef_XT(ixt,i)
679             write(*,*) 'Bcoef_XT(ixt,i)=',Bcoef_XT(ixt,i)
680             write(*,*) 'flx_xt1(ixt,i)=',flx_xt1(ixt,i)
681             if (k.ge.2) then
682               write(*,*) 'Ccoef_XT(ixt,i,k)=',Ccoef_XT(ixt,i,k)
683               write(*,*) 'Dcoef_XT(ixt,i,k)=',Dcoef_XT(ixt,i,k)
684             endif
685             stop
686         endif
687        enddo !do ixt=1,ntraciso
688       END DO
689    END DO       
690#endif
691#ifdef ISOVERIF
692     if (iso_eau.gt.0) then
693        call iso_verif_egalite_vect2D( &
694     &           xt_new,q_new, &
695     &           'climb_hq_mod 504',ntraciso,klon,klev)
696     endif !if (iso_eau.gt.0) then
697     if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then 
698        do i=1,klon
699         do k=1,klev
700            if (q_new(i,k).gt.ridicule) then 
701               if (iso_verif_o18_aberrant_nostop( &
702     &              xt_new(iso_HDO,i,k)/q_new(i,k), &
703     &              xt_new(iso_O18,i,k)/q_new(i,k), &
704     &              'climb_hq_mod 690').eq.1) then
705                  write(*,*) 'i,k,q_new(i,k)=',i,k,q_new(i,k)                       
706                  stop
707              endif !  if (iso_verif_o18_aberrant_nostop
708            endif !if (q_seri(i,k).gt.errmax) then   
709          enddo !k=1,klev
710         enddo  !i=1,klon
711        endif !if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then   
712#endif
713#endif   
714!****************************************************************************************
715! 3)
716! Calculation of the flux for Q and H
717!
718!****************************************************************************************
719
720!- The flux at first layer, k=1
721    flux_q(1:knon,1)=flx_q1(1:knon)
722    flux_h(1:knon,1)=flx_h1(1:knon)
723#ifdef ISO
724    do ixt=1,ntraciso
725        flux_xt(ixt,1:knon,1)=flx_xt1(ixt,1:knon)
726    enddo ! do ixt=1,ntraciso
727#endif
728
729!- The flux at all layers above surface
730    DO k = 2, klev
731       DO i = 1, knon
732          flux_q(i,k) = (Kcoefhq(i,k)/RG/dtime) * &
733               (q_new(i,k)-q_new(i,k-1)+gamaq(i,k))
734
735          flux_h(i,k) = (Kcoefhq(i,k)/RG/dtime) * &
736               (h_new(i,k)-h_new(i,k-1)+gamah(i,k))
737
738#ifdef ISO
739        do ixt=1,ntraciso
740          flux_xt(ixt,i,k) = (Kcoefhq(i,k)/RG/dtime) * &
741               (xt_new(ixt,i,k)-xt_new(ixt,i,k-1)+gamaxt(ixt,i,k))
742        enddo ! do ixt=1,ntraciso
743#endif
744       END DO
745    END DO
746
747!****************************************************************************************
748! 4)
749! Calculation of tendency for Q and H
750!
751!****************************************************************************************
752    d_h_col_vdf(:) = 0.0
753    DO k = 1, klev
754       DO i = 1, knon
755          d_t(i,k) = h_new(i,k)/(psref(i)/pplay(i,k))**RKAPPA/RCPD - t_old(i,k)
756          d_q(i,k) = q_new(i,k) - q_old(i,k)
757          d_h(i,k) = h_new(i,k) - h_old(i,k)
758!JLD          d_t(i,k) = d_h(i,k)/(psref(i)/pplay(i,k))**RKAPPA/RCPD    !correction a venir
759    ! layer air mass
760          zairm(i, k) = (paprs(i,k)-paprs(i,k+1))/rg
761          d_h_col_vdf(i) = d_h_col_vdf(i) + d_h(i,k)*zairm(i,k)
762#ifdef ISO
763        do ixt=1,ntraciso
764          d_xt(ixt,i,k) = xt_new(ixt,i,k) - xt_old(ixt,i,k)
765        enddo ! do ixt=1,ntraciso
766#ifdef ISOVERIF
767        do ixt=1,ntraciso
768          call iso_verif_noNaN(xt_new(ixt,i,k),'climb_hq 562')
769          call iso_verif_noNaN(d_xt(ixt,i,k),'climb_hq 562')
770        enddo ! do ixt=1,ntraciso
771#endif
772#ifdef ISOVERIF
773        if (iso_eau.gt.0) then
774           call iso_verif_egalite(d_xt(iso_eau,i,k), &
775     &         d_q(i,k),'climb_hq_mod 503')
776           call iso_verif_egalite(xt_new(iso_eau,i,k), &
777     &         q_new(i,k),'climb_hq_mod 503b')
778        endif
779#ifdef ISOTRAC
780        IF(nzone > 0) CALL iso_verif_traceur(xt_old(1,i,k),'climb_hq_mod 526')
781#endif       
782#endif       
783#endif
784        END DO
785    END DO
786
787#ifdef ISO
788#ifdef ISOVERIF
789!        write(*,*) 'climb_hq_mod 758: d_xt,d_q=',d_xt(iso_eau,1,1),d_q(1,1)
790        if (iso_eau.gt.0) then
791         call iso_verif_egalite_vect2D( &
792                d_xt,d_q, &
793                'climb_hq_mod 761',ntraciso,klon,klev)
794        endif       
795#endif
796#endif
797
798!****************************************************************************************
799! Some deallocations
800!
801!****************************************************************************************
802    IF (last) THEN
803       DEALLOCATE(Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H,stat=ierr)   
804       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H, ierr=', ierr
805       DEALLOCATE(Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H,stat=ierr)   
806       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H, ierr=', ierr
807       DEALLOCATE(gamaq, gamah,stat=ierr)
808       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate gamaq, gamah, ierr=', ierr
809       DEALLOCATE(Kcoefhq,stat=ierr)
810       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Kcoefhq, ierr=', ierr
811       DEALLOCATE(h_old, d_h_col_vdf, f_h_bnd, stat=ierr)
812       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate h_old, d_h_col_vdf, f_h_bnd, ierr=', ierr
813
814#ifdef ISO
815       DEALLOCATE(Ccoef_XT, Dcoef_XT ,stat=ierr)   
816       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Ccoef_XT, Dcoef_XT, ierr=', ierr
817       DEALLOCATE(Acoef_XT, Bcoef_XT, stat=ierr)   
818       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Acoef_XT, Bcoef_XT, ierr=', ierr
819       DEALLOCATE(gamaxt, stat=ierr)
820       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate gamaxt, ierr=', ierr
821#endif
822    END IF
823  END SUBROUTINE climb_hq_up
824!
825!****************************************************************************************
826!
827END MODULE climb_hq_mod
828
829 
830
831
832
833
Note: See TracBrowser for help on using the repository browser.