source: LMDZ6/trunk/libf/phylmdiso/ajsec.F90 @ 5894

Last change on this file since 5894 was 5894, checked in by Sebastien Nguyen, 3 weeks ago

rephase LMDZISO with 5864 version of phylmd + bug fixes in physiq_mod + other bugs in isoverif sections. Code now compiles and runs with -debug -isotopes true -isoverif. There are still isoverif error messages for Dexcess getting greater than 1000 on some points at some moments.
  • Property svn:keywords set to Id
File size: 15.3 KB
Line 
1
2! $Header$
3!$gpum horizontal klon
4MODULE ajsec_mod
5  PRIVATE
6
7  PUBLIC ajsec, ajsec_convv2, ajsec_old
8
9  CONTAINS
10
11SUBROUTINE ajsec(paprs, pplay, t, q, limbas, d_t, d_q &
12#ifdef ISO     
13               ,xt,d_xt &
14#endif
15          )
16  USE dimphy
17#ifdef ISO
18  USE infotrac_phy, ONLY: ntraciso =>ntiso   
19#ifdef ISOVERIF
20  USE isotopes_mod, ONLY : iso_eau,iso_HDO
21  USE isotopes_verif_mod, ONLY: iso_verif_egalite, &
22        iso_verif_egalite_choix,iso_verif_noNaN,errmax,errmaxrel
23#ifdef ISOTRAC
24  USE isotopes_verif_mod, ONLY: iso_verif_traceur,iso_verif_traceur_justmass
25#endif
26#endif
27#endif
28  USE yomcst_mod_h
29IMPLICIT NONE
30  ! ======================================================================
31  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
32  ! Objet: ajustement sec (adaptation du GCM du LMD)
33  ! ======================================================================
34  ! Arguments:
35  ! t-------input-R- Temperature
36
37  ! d_t-----output-R-Incrementation de la temperature
38  ! ======================================================================
39
40  REAL paprs(klon, klev+1), pplay(klon, klev)
41  REAL t(klon, klev), q(klon, klev)
42  REAL d_t(klon, klev), d_q(klon, klev)
43
44  INTEGER limbas(klon), limhau ! les couches a ajuster
45
46  LOGICAL mixq
47  ! cc      PARAMETER (mixq=.TRUE.)
48  PARAMETER (mixq=.FALSE.)
49
50  REAL zh(klon, klev)
51  REAL zho(klon, klev)
52  REAL zq(klon, klev)
53  REAL zpk(klon, klev)
54  REAL zpkdp(klon, klev)
55  REAL hm, sm, qm
56  LOGICAL modif(klon), down
57  INTEGER i, k, k1, k2
58#ifdef ISO
59  REAL xt(ntraciso,klon,klev)
60  REAL d_xt(ntraciso,klon,klev)
61  REAL zxt(ntraciso,klon,klev)
62  REAL xtm(ntraciso)
63  INTEGER ixt
64#endif
65
66  ! Initialisation:
67
68#ifdef ISO
69#ifdef ISOVERIF
70  DO i=1,klon
71    DO k=1,klev         
72      IF (iso_eau.GT.0) THEN
73        CALL iso_verif_egalite_choix(q(i,k),xt(iso_eau,i,k), &
74             'ajsec 67',errmax,errmaxrel)
75      ENDIF !if (iso_eau.gt.0) then
76    ENDDO !do k=limbas,limhau
77  ENDDO !do i=1,klon
78#endif
79#endif
80
81  ! ym
82  limhau = klev
83
84  DO k = 1, klev
85    DO i = 1, klon
86      d_t(i, k) = 0.0
87      d_q(i, k) = 0.0
88#ifdef ISO
89      DO ixt=1,ntraciso
90         d_xt(ixt,i,k)=0.0
91      ENDDO
92#endif
93    ENDDO
94  ENDDO
95  ! ------------------------------------- detection des profils a modifier
96  DO k = 1, limhau
97    DO i = 1, klon
98      zpk(i, k) = pplay(i, k)**rkappa
99      zh(i, k) = rcpd*t(i, k)/zpk(i, k)
100      zho(i, k) = zh(i, k)
101      zq(i, k) = q(i, k)
102#ifdef ISO
103      DO ixt=1,ntraciso
104         zxt(ixt,i,k)=xt(ixt,i,k)
105      ENDDO
106#endif
107    ENDDO
108  ENDDO
109
110  DO k = 1, limhau
111    DO i = 1, klon
112      zpkdp(i, k) = zpk(i, k)*(paprs(i,k)-paprs(i,k+1))
113    ENDDO
114  ENDDO
115
116  DO i = 1, klon
117    modif(i) = .FALSE.
118  ENDDO
119  DO k = 2, limhau
120    DO i = 1, klon
121      IF (.NOT. modif(i) .AND. k-1>limbas(i)) THEN
122        IF (zh(i,k)<zh(i,k-1)) modif(i) = .TRUE.
123      ENDIF
124    ENDDO
125  ENDDO
126  ! ------------------------------------- correction des profils instables
127  DO i = 1, klon
128    IF (modif(i)) THEN
129      k2 = limbas(i)
1308000  CONTINUE
131      k2 = k2 + 1
132      IF (k2>limhau) GO TO 8001
133      IF (zh(i,k2)<zh(i,k2-1)) THEN
134        k1 = k2 - 1
135        k = k1
136        sm = zpkdp(i, k2)
137        hm = zh(i, k2)
138        qm = zq(i, k2)
139#ifdef ISO
140        DO ixt=1,ntraciso
141           xtm(ixt)=zxt(ixt,i,k2)
142        ENDDO   
143#ifdef ISOVERIF
144        IF (iso_eau.GT.0) THEN
145           CALL iso_verif_egalite_choix(qm,xtm(iso_eau), &
146                             'ajsec 126',errmax,errmaxrel)
147        ENDIF !if (iso_eau.gt.0) then     
148#endif
149#endif
1508020    CONTINUE
151        sm = sm + zpkdp(i, k)
152        hm = hm + zpkdp(i, k)*(zh(i,k)-hm)/sm
153        qm = qm + zpkdp(i, k)*(zq(i,k)-qm)/sm
154#ifdef ISO
155        DO ixt=1,ntraciso
156           xtm(ixt)=xtm(ixt) &
157                   +zpkdp(i,k)*(zxt(ixt,i,k)-xtm(ixt))/sm
158        ENDDO   
159#ifdef ISOVERIF
160        IF (iso_eau.GT.0) THEN
161           CALL iso_verif_egalite_choix(qm,xtm(iso_eau), &
162                             'ajsec 136',errmax,errmaxrel)
163        ENDIF !if (iso_eau.gt.0) then     
164#endif
165#endif
166        down = .FALSE.
167        IF (k1/=limbas(i)) THEN
168          IF (hm<zh(i,k1-1)) down = .TRUE.
169        ENDIF
170        IF (down) THEN
171          k1 = k1 - 1
172          k = k1
173        ELSE
174          IF ((k2==limhau)) GO TO 8021
175          IF ((zh(i,k2+1)>=hm)) GO TO 8021
176          k2 = k2 + 1
177          k = k2
178        ENDIF
179        GO TO 8020
1808021    CONTINUE
181        ! ------------ nouveau profil : constant (valeur moyenne)
182        DO k = k1, k2
183          zh(i, k) = hm
184          zq(i, k) = qm
185#ifdef ISO
186          DO ixt=1,ntraciso
187             zxt(ixt,i,k)=xtm(ixt)
188          ENDDO
189#endif
190        ENDDO
191        k2 = k2 + 1
192      ENDIF
193      GO TO 8000
1948001  CONTINUE
195    ENDIF
196  ENDDO
197
198#ifdef ISO
199      ! cam verif
200#ifdef ISOVERIF
201  DO i=1,klon
202    DO k=1,klev
203      DO ixt=1,ntraciso
204        CALL iso_verif_noNaN(zxt(ixt,i,k),'ajsec 173')
205      ENDDO !do ixt=1,niso           
206      IF (iso_eau.GT.0) THEN
207        CALL iso_verif_egalite_choix(zq(i,k),zxt(iso_eau,i,k), &
208             'ajsec 168',errmax,errmaxrel)
209      ENDIF !if (iso_eau.gt.0) then
210#ifdef ISOTRAC     
211      CALL iso_verif_traceur(zxt(1,i,k),'ajsec 181')
212#endif           
213    ENDDO !do k=limbas,limhau
214  ENDDO !do i=1,klon
215#endif
216      ! end cam verif
217#endif
218
219  DO k = 1, limhau
220    DO i = 1, klon
221      d_t(i, k) = (zh(i,k)-zho(i,k))*zpk(i, k)/rcpd
222      d_q(i, k) = zq(i, k) - q(i, k)
223#ifdef ISO
224      DO ixt=1,ntraciso
225        d_xt(ixt,i,k)=zxt(ixt,i,k)-xt(ixt,i,k)
226      ENDDO
227#endif
228    ENDDO
229  ENDDO
230
231#ifdef ISO
232      ! cam verif
233#ifdef ISOVERIF
234  DO i = 1, klon
235    DO k = 1, limhau
236      IF (iso_eau.GT.0) THEN
237        CALL iso_verif_egalite_choix(d_q(i,k),d_xt(iso_eau,i,k), &
238             'ajsec 198',errmax,errmaxrel)
239      ENDIF
240#ifdef ISOTRAC     
241      CALL iso_verif_traceur_justmass(d_xt(1,i,k),'physiq 210')
242#endif         
243    ENDDO
244  ENDDO
245#endif
246      ! end cam verif     
247#endif
248
249  ! FH : les d_q et d_t sont maintenant calcules de facon a valoir
250  ! effectivement 0. si on ne fait rien.
251
252  ! IF (limbas.GT.1) THEN
253  ! DO k = 1, limbas-1
254  ! DO i = 1, klon
255  ! d_t(i,k) = 0.0
256  ! d_q(i,k) = 0.0
257  ! ENDDO
258  ! ENDDO
259  ! ENDIF
260
261  ! IF (limhau.LT.klev) THEN
262  ! DO k = limhau+1, klev
263  ! DO i = 1, klon
264  ! d_t(i,k) = 0.0
265  ! d_q(i,k) = 0.0
266  ! ENDDO
267  ! ENDDO
268  ! ENDIF
269
270  IF (.NOT. mixq) THEN
271    DO k = 1, klev
272      DO i = 1, klon
273        d_q(i, k) = 0.0
274#ifdef ISO
275        DO ixt=1,ntraciso
276           d_xt(ixt,i,k)=0.0
277        ENDDO
278#endif
279      ENDDO
280    ENDDO
281  ENDIF
282
283#ifdef ISO
284      ! cam verif
285#ifdef ISOVERIF
286  DO i = 1, klon
287    DO k = 1, klev
288      IF (iso_eau.GT.0) THEN
289        CALL iso_verif_egalite(d_q(i,k),d_xt(iso_eau,i,k),'ajsec 270')
290      ENDIF
291#ifdef ISOTRAC     
292      CALL iso_verif_traceur_justmass(d_xt(1,i,k),'physiq 3045')
293#endif         
294    ENDDO
295  ENDDO
296#endif
297      ! end cam verif
298#endif
299
300  RETURN
301END SUBROUTINE ajsec
302
303SUBROUTINE ajsec_convv2(paprs, pplay, t, q, d_t, d_q &
304#ifdef ISO
305           ,xt,d_xt &     
306#endif         
307        )
308  USE dimphy
309#ifdef ISO
310  USE infotrac_phy, ONLY: ntraciso=>ntiso   
311#ifdef ISOVERIF
312  USE isotopes_mod, ONLY : iso_eau,iso_HDO
313  USE isotopes_verif_mod, ONLY: iso_verif_egalite, &
314        iso_verif_egalite_choix,iso_verif_noNaN,errmax,errmaxrel
315#ifdef ISOTRAC
316  USE isotopes_verif_mod, ONLY: iso_verif_traceur,iso_verif_traceur_justmass
317#endif
318#endif
319#endif
320  USE yomcst_mod_h
321IMPLICIT NONE
322  ! ======================================================================
323  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
324  ! Objet: ajustement sec (adaptation du GCM du LMD)
325  ! ======================================================================
326  ! Arguments:
327  ! t-------input-R- Temperature
328
329  ! d_t-----output-R-Incrementation de la temperature
330  ! ======================================================================
331
332  REAL paprs(klon, klev+1), pplay(klon, klev)
333  REAL t(klon, klev), q(klon, klev)
334  REAL d_t(klon, klev), d_q(klon, klev)
335
336  INTEGER limbas, limhau ! les couches a ajuster
337  ! cc      PARAMETER (limbas=klev-3, limhau=klev)
338  ! ym      PARAMETER (limbas=1, limhau=klev)
339
340  LOGICAL mixq
341  ! cc      PARAMETER (mixq=.TRUE.)
342  PARAMETER (mixq=.FALSE.)
343
344  REAL zh(klon, klev)
345  REAL zq(klon, klev)
346  REAL zpk(klon, klev)
347  REAL zpkdp(klon, klev)
348  REAL hm, sm, qm
349  LOGICAL modif(klon), down
350  INTEGER i, k, k1, k2
351
352#ifdef ISO
353  REAL xt(ntraciso,klon,klev)
354  REAL d_xt(ntraciso,klon,klev)
355  REAL zxt(ntraciso,klon,klev)
356  REAL xtm(ntraciso)
357  INTEGER ixt
358#endif
359
360
361#ifdef ISO
362  ! cam verif
363#ifdef ISOVERIF
364  DO i=1,klon
365    DO k=1,klev
366      DO ixt=1,ntraciso
367        CALL iso_verif_noNAN(xt(ixt,i,k),'ajsec 320')
368      ENDDO !do ixt=1,niso           
369      IF (iso_eau.GT.0) THEN
370        CALL iso_verif_egalite_choix(q(i,k),xt(iso_eau,i,k), &
371             'ajsec 324',errmax,errmaxrel)
372      ENDIF !if (iso_eau.gt.0) then
373#ifdef ISOTRAC     
374      CALL iso_verif_traceur(xt(1,i,k),'ajsec 327')
375#endif           
376    ENDDO !do k=1,klev
377  ENDDO !do i=1,klon
378#endif
379      ! end cam verif
380#endif
381
382  ! Initialisation:
383
384  ! ym
385  limbas = 1
386  limhau = klev
387
388  DO k = 1, klev
389    DO i = 1, klon
390      d_t(i, k) = 0.0
391      d_q(i, k) = 0.0
392#ifdef ISO
393      DO ixt=1,ntraciso
394        d_xt(ixt,i,k)=0.0
395      ENDDO
396#endif
397    ENDDO
398  ENDDO
399  ! ------------------------------------- detection des profils a modifier
400  DO k = limbas, limhau
401    DO i = 1, klon
402      zpk(i, k) = pplay(i, k)**rkappa
403      zh(i, k) = rcpd*t(i, k)/zpk(i, k)
404      zq(i, k) = q(i, k)
405#ifdef ISO
406      DO ixt=1,ntraciso
407        zxt(ixt,i,k)=xt(ixt,i,k)
408      ENDDO
409#endif
410    ENDDO
411  ENDDO
412
413  DO k = limbas, limhau
414    DO i = 1, klon
415      zpkdp(i, k) = zpk(i, k)*(paprs(i,k)-paprs(i,k+1))
416    ENDDO
417  ENDDO
418
419  DO i = 1, klon
420    modif(i) = .FALSE.
421  ENDDO
422  DO k = limbas + 1, limhau
423    DO i = 1, klon
424      IF (.NOT. modif(i)) THEN
425        IF (zh(i,k)<zh(i,k-1)) modif(i) = .TRUE.
426      ENDIF
427    ENDDO
428  ENDDO
429  ! ------------------------------------- correction des profils instables
430  DO i = 1, klon
431    IF (modif(i)) THEN
432      k2 = limbas
4338000  CONTINUE
434      k2 = k2 + 1
435      IF (k2>limhau) GO TO 8001
436      IF (zh(i,k2)<zh(i,k2-1)) THEN
437        k1 = k2 - 1
438        k = k1
439        sm = zpkdp(i, k2)
440        hm = zh(i, k2)
441        qm = zq(i, k2)
442#ifdef ISO
443        DO ixt=1,ntraciso
444          xtm(ixt)=zxt(ixt,i,k2)
445        ENDDO
446#endif
4478020    CONTINUE
448        sm = sm + zpkdp(i, k)
449        hm = hm + zpkdp(i, k)*(zh(i,k)-hm)/sm
450        qm = qm + zpkdp(i, k)*(zq(i,k)-qm)/sm
451#ifdef ISO
452        DO ixt=1,ntraciso
453          xtm(ixt)=xtm(ixt) &
454                  +zpkdp(i,k)*(zxt(ixt,i,k)-xtm(ixt))/sm
455        ENDDO
456#endif
457        down = .FALSE.
458        IF (k1/=limbas) THEN
459          IF (hm<zh(i,k1-1)) down = .TRUE.
460        ENDIF
461        IF (down) THEN
462          k1 = k1 - 1
463          k = k1
464        ELSE
465          IF ((k2==limhau)) GO TO 8021
466          IF ((zh(i,k2+1)>=hm)) GO TO 8021
467          k2 = k2 + 1
468          k = k2
469        ENDIF
470        GO TO 8020
4718021    CONTINUE
472        ! ------------ nouveau profil : constant (valeur moyenne)
473        DO k = k1, k2
474          zh(i, k) = hm
475          zq(i, k) = qm
476#ifdef ISO
477          DO ixt=1,ntraciso
478            zxt(ixt,i,k)=xtm(ixt)
479          ENDDO
480#endif
481        ENDDO
482        k2 = k2 + 1
483      ENDIF
484      GO TO 8000
4858001  CONTINUE
486    ENDIF
487  ENDDO
488
489#ifdef ISO
490      ! cam verif
491#ifdef ISOVERIF
492  DO i=1,klon
493    DO k=limbas,limhau
494      DO ixt=1,ntraciso
495        CALL iso_verif_noNAN(zxt(ixt,i,k),'ajsec 428')
496      ENDDO !do ixt=1,niso           
497      IF (iso_eau.GT.0) THEN
498        CALL iso_verif_egalite_choix(zq(i,k),zxt(iso_eau,i,k), &
499             'ajsec 432',errmax,errmaxrel)
500      ENDIF !if (iso_eau.gt.0) then
501#ifdef ISOTRAC     
502      CALL iso_verif_traceur(zxt(1,i,k),'ajsec 436')
503#endif           
504    ENDDO !do k=limbas,limhau
505  ENDDO !do i=1,klon
506#endif
507      ! end cam verif
508#endif   
509
510
511  DO k = limbas, limhau
512    DO i = 1, klon
513      d_t(i, k) = zh(i, k)*zpk(i, k)/rcpd - t(i, k)
514      d_q(i, k) = zq(i, k) - q(i, k)
515#ifdef ISO
516      DO ixt=1,ntraciso
517         d_xt(ixt,i,k)=zxt(ixt,i,k)-xt(ixt,i,k)
518      ENDDO
519#endif
520    ENDDO
521  ENDDO
522
523  IF (limbas>1) THEN
524    DO k = 1, limbas - 1
525      DO i = 1, klon
526        d_t(i, k) = 0.0
527        d_q(i, k) = 0.0
528#ifdef ISO
529        DO ixt=1,ntraciso
530          d_xt(ixt,i,k)=0.0
531        ENDDO
532#endif
533      ENDDO
534    ENDDO
535  ENDIF
536
537  IF (limhau<klev) THEN
538    DO k = limhau + 1, klev
539      DO i = 1, klon
540        d_t(i, k) = 0.0
541        d_q(i, k) = 0.0
542#ifdef ISO
543        DO ixt=1,ntraciso
544          d_xt(ixt,i,k)=0.0
545        ENDDO
546#endif
547      ENDDO
548    ENDDO
549  ENDIF
550
551  IF (.NOT. mixq) THEN
552    DO k = 1, klev
553      DO i = 1, klon
554        d_q(i, k) = 0.0
555#ifdef ISO
556        DO ixt=1,ntraciso
557          d_xt(ixt,i,k)=0.0
558        ENDDO
559#endif
560      ENDDO
561    ENDDO
562  ENDIF
563
564
565#ifdef ISO
566      ! cam verif
567#ifdef ISOVERIF
568  DO i = 1, klon
569    DO k = limbas, limhau
570      IF (iso_eau.GT.0) THEN
571        CALL iso_verif_egalite(d_q(i,k),d_xt(iso_eau,i,k),'ajsec 270')
572      ENDIF
573#ifdef ISOTRAC     
574      CALL iso_verif_traceur_justmass(d_xt(1,i,k),'physiq 3045')
575#endif         
576    ENDDO
577  ENDDO
578#endif
579      ! end cam verif
580#endif
581
582  RETURN
583END SUBROUTINE ajsec_convv2
584
585SUBROUTINE ajsec_old(paprs, pplay, t, d_t)
586  USE dimphy
587  USE yomcst_mod_h
588  IMPLICIT NONE
589  ! ======================================================================
590  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
591  ! Objet: ajustement sec (adaptation du GCM du LMD)
592  ! ======================================================================
593  ! Arguments:
594  ! t-------input-R- Temperature
595
596  ! d_t-----output-R-Incrementation de la temperature
597  ! ======================================================================
598
599  REAL paprs(klon, klev+1), pplay(klon, klev)
600  REAL t(klon, klev)
601  REAL d_t(klon, klev)
602
603  REAL local_h(klon, klev)
604  REAL hm, sm
605  LOGICAL modif(klon), down
606  INTEGER i, l, l1, l2
607  ! ------------------------------------- detection des profils a modifier
608  DO i = 1, klon
609    modif(i) = .FALSE.
610  END DO
611
612  DO l = 1, klev
613    DO i = 1, klon
614      local_h(i, l) = rcpd*t(i, l)/(pplay(i,l)**rkappa)
615    END DO
616  END DO
617
618  DO l = 2, klev
619    DO i = 1, klon
620      IF (local_h(i,l)<local_h(i,l-1)) THEN
621        modif(i) = .TRUE.
622      ELSE
623        modif(i) = modif(i)
624      END IF
625    END DO
626  END DO
627  ! ------------------------------------- correction des profils instables
628  DO i = 1, klon
629    IF (modif(i)) THEN
630      l2 = 1
6318000  CONTINUE
632      l2 = l2 + 1
633      IF (l2>klev) GO TO 8001
634      IF (local_h(i,l2)<local_h(i,l2-1)) THEN
635        l1 = l2 - 1
636        l = l1
637        sm = pplay(i, l2)**rkappa*(paprs(i,l2)-paprs(i,l2+1))
638        hm = local_h(i, l2)
6398020    CONTINUE
640        sm = sm + pplay(i, l)**rkappa*(paprs(i,l)-paprs(i,l+1))
641        hm = hm + pplay(i, l)**rkappa*(paprs(i,l)-paprs(i,l+1))*(local_h(i,l) &
642          -hm)/sm
643        down = .FALSE.
644        IF (l1/=1) THEN
645          IF (hm<local_h(i,l1-1)) THEN
646            down = .TRUE.
647          END IF
648        END IF
649        IF (down) THEN
650          l1 = l1 - 1
651          l = l1
652        ELSE
653          IF ((l2==klev)) GO TO 8021
654          IF ((local_h(i,l2+1)>=hm)) GO TO 8021
655          l2 = l2 + 1
656          l = l2
657        END IF
658        GO TO 8020
6598021    CONTINUE
660        ! ------------ nouveau profil : constant (valeur moyenne)
661        DO l = l1, l2
662          local_h(i, l) = hm
663        END DO
664        l2 = l2 + 1
665      END IF
666      GO TO 8000
6678001  CONTINUE
668    END IF
669  END DO
670
671  DO l = 1, klev
672    DO i = 1, klon
673      d_t(i, l) = local_h(i, l)*(pplay(i,l)**rkappa)/rcpd - t(i, l)
674    END DO
675  END DO
676
677  RETURN
678END SUBROUTINE ajsec_old
679
680END MODULE ajsec_mod
Note: See TracBrowser for help on using the repository browser.