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

Last change on this file since 5837 was 5285, checked in by abarral, 13 months ago

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