source: LMDZ6/branches/Amaury_dev/libf/phylmdiso/ajsec.F90 @ 5104

Last change on this file since 5104 was 5103, checked in by abarral, 4 months ago

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