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

Last change on this file since 5213 was 5158, checked in by abarral, 3 months ago

Add missing klon on strataer_emiss_mod.F90
Correct various missing explicit declarations
Replace tabs by spaces (tabs are not part of the fortran charset)
Continue cleaning modules
Removed unused arguments and variables
  • Property svn:keywords set to Id
File size: 15.6 KB
RevLine 
[3927]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
[4143]11    USE infotrac_phy, ONLY: ntraciso =>ntiso   
[3927]12#ifdef ISOVERIF
[5101]13  USE isotopes_mod, ONLY: iso_eau,iso_HDO
[3927]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
[5144]21USE lmdz_yomcst
[3927]22  IMPLICIT 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  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
[5117]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
[3927]56#endif
57
58
59  ! Initialisation:
60
61#ifdef ISO
62#ifdef ISOVERIF
[5158]63      DO i=1,klon
64         DO k=1,klev
[5117]65           IF (iso_eau.gt.0) THEN
[5103]66             CALL iso_verif_egalite_choix(q(i,k),xt(iso_eau,i,k), &
[3927]67                 'ajsec 67',errmax,errmaxrel)
[5116]68           endif !if (iso_eau.gt.0) THEN
[3927]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
[5158]83         DO ixt=1,ntraciso
[3927]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
[5158]97      DO ixt=1,ntraciso
[3927]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
[5158]134        DO ixt=1,ntraciso
[3927]135                 xtm(ixt)=zxt(ixt,i,k2)
136        enddo   
137#ifdef ISOVERIF
[5117]138        IF (iso_eau.gt.0) THEN
[5103]139             CALL iso_verif_egalite_choix(qm,xtm(iso_eau), &
[3927]140                 'ajsec 126',errmax,errmaxrel)
[5116]141        endif !if (iso_eau.gt.0) THEN
[3927]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
[5158]149        DO ixt=1,ntraciso
[3927]150                 xtm(ixt)=xtm(ixt) &
151                    +zpkdp(i,k)*(zxt(ixt,i,k)-xtm(ixt))/sm
152        enddo   
153#ifdef ISOVERIF
[5117]154        IF (iso_eau.gt.0) THEN
[5103]155             CALL iso_verif_egalite_choix(qm,xtm(iso_eau), &
[3927]156                 'ajsec 136',errmax,errmaxrel)
[5116]157        endif !if (iso_eau.gt.0) THEN
[3927]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
[5158]180                DO ixt=1,ntraciso
[3927]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
[5158]195      DO i=1,klon
196         DO k=1,klev
197           DO ixt=1,ntraciso
[5103]198             CALL iso_verif_noNaN(zxt(ixt,i,k),'ajsec 173')
[3927]199           enddo !do ixt=1,niso           
[5117]200           IF (iso_eau.gt.0) THEN
[5103]201             CALL iso_verif_egalite_choix(zq(i,k),zxt(iso_eau,i,k), &
[3927]202                 'ajsec 168',errmax,errmaxrel)
[5116]203           endif !if (iso_eau.gt.0) THEN
[3927]204#ifdef ISOTRAC     
[5103]205           CALL iso_verif_traceur(zxt(1,i,k),'ajsec 181')
[3927]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
[5158]218         DO ixt=1,ntraciso
[3927]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
[5158]228      DO i = 1, klon
229        DO k = 1, limhau
[5117]230         IF (iso_eau.gt.0) THEN
[5103]231          CALL iso_verif_egalite_choix(d_q(i,k),d_xt(iso_eau,i,k), &
[3927]232                'ajsec 198',errmax,errmaxrel)
233         endif
234#ifdef ISOTRAC     
[5103]235        CALL iso_verif_traceur_justmass(d_xt(1,i,k),'physiq 210')
[3927]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
[5158]269         DO ixt=1,ntraciso
[3927]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
[5158]280      DO i = 1, klon
281        DO k = 1, klev
[5117]282         IF (iso_eau.gt.0) THEN
[5103]283          CALL iso_verif_egalite(d_q(i,k),d_xt(iso_eau,i,k),'ajsec 270')
[3927]284         endif
285#ifdef ISOTRAC     
[5103]286        CALL iso_verif_traceur_justmass(d_xt(1,i,k),'physiq 3045')
[3927]287#endif         
288         enddo
289      enddo
290#endif
291      ! end cam verif
292#endif
293
294
[5105]295
[3927]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
[4143]305    USE infotrac_phy, ONLY: ntraciso=>ntiso   
[3927]306#ifdef ISOVERIF
[5101]307  USE isotopes_mod, ONLY: iso_eau,iso_HDO
[3927]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
[5144]315        USE lmdz_yomcst
316
[3927]317  IMPLICIT 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  REAL paprs(klon, klev+1), pplay(klon, klev)
328  REAL t(klon, klev), q(klon, klev)
329  REAL d_t(klon, klev), d_q(klon, klev)
330
331  INTEGER limbas, limhau ! les couches a ajuster
332  ! cc      PARAMETER (limbas=klev-3, limhau=klev)
333  ! ym      PARAMETER (limbas=1, limhau=klev)
334
335  LOGICAL mixq
336  ! cc      PARAMETER (mixq=.TRUE.)
337  PARAMETER (mixq=.FALSE.)
338
339  REAL zh(klon, klev)
340  REAL zq(klon, klev)
341  REAL zpk(klon, klev)
342  REAL zpkdp(klon, klev)
343  REAL hm, sm, qm
344  LOGICAL modif(klon), down
345  INTEGER i, k, k1, k2
346
347#ifdef ISO
[5117]348      REAL xt(ntraciso,klon,klev)
349      REAL d_xt(ntraciso,klon,klev)
350      REAL zxt(ntraciso,klon,klev)
351      REAL xtm(ntraciso)
352      INTEGER ixt
[3927]353#endif
354
355
356#ifdef ISO
357      ! cam verif
358#ifdef ISOVERIF
[5158]359      DO i=1,klon
360         DO k=1,klev
361           DO ixt=1,ntraciso
[5103]362             CALL iso_verif_noNAN(xt(ixt,i,k),'ajsec 320')
[3927]363           enddo !do ixt=1,niso           
[5117]364           IF (iso_eau.gt.0) THEN
[5103]365             CALL iso_verif_egalite_choix(q(i,k),xt(iso_eau,i,k), &
[3927]366                 'ajsec 324',errmax,errmaxrel)
[5116]367           endif !if (iso_eau.gt.0) THEN
[3927]368#ifdef ISOTRAC     
[5103]369           CALL iso_verif_traceur(xt(1,i,k),'ajsec 327')
[3927]370#endif           
371         enddo !do k=1,klev
372      enddo !do i=1,klon
373#endif
374      ! end cam verif
375#endif
376
377  ! Initialisation:
378
379  ! ym
380  limbas = 1
381  limhau = klev
382
383  DO k = 1, klev
384    DO i = 1, klon
385      d_t(i, k) = 0.0
386      d_q(i, k) = 0.0
387#ifdef ISO
[5158]388         DO ixt=1,ntraciso
[3927]389            d_xt(ixt,i,k)=0.0
390         enddo
391#endif
392    END DO
393  END DO
394  ! ------------------------------------- detection des profils a modifier
395  DO k = limbas, limhau
396    DO i = 1, klon
397      zpk(i, k) = pplay(i, k)**rkappa
398      zh(i, k) = rcpd*t(i, k)/zpk(i, k)
399      zq(i, k) = q(i, k)
400#ifdef ISO
[5158]401         DO ixt=1,ntraciso
[3927]402            zxt(ixt,i,k)=xt(ixt,i,k)
403         enddo
404#endif
405
406    END DO
407  END DO
408
409  DO k = limbas, limhau
410    DO i = 1, klon
411      zpkdp(i, k) = zpk(i, k)*(paprs(i,k)-paprs(i,k+1))
412    END DO
413  END DO
414
415  DO i = 1, klon
416    modif(i) = .FALSE.
417  END DO
418  DO k = limbas + 1, limhau
419    DO i = 1, klon
420      IF (.NOT. modif(i)) THEN
421        IF (zh(i,k)<zh(i,k-1)) modif(i) = .TRUE.
422      END IF
423    END DO
424  END DO
425  ! ------------------------------------- correction des profils instables
426  DO i = 1, klon
427    IF (modif(i)) THEN
428      k2 = limbas
4298000  CONTINUE
430      k2 = k2 + 1
431      IF (k2>limhau) GO TO 8001
432      IF (zh(i,k2)<zh(i,k2-1)) THEN
433        k1 = k2 - 1
434        k = k1
435        sm = zpkdp(i, k2)
436        hm = zh(i, k2)
437        qm = zq(i, k2)
438#ifdef ISO
[5158]439              DO ixt=1,ntraciso
[3927]440                 xtm(ixt)=zxt(ixt,i,k2)
441              enddo
442#endif
4438020    CONTINUE
444        sm = sm + zpkdp(i, k)
445        hm = hm + zpkdp(i, k)*(zh(i,k)-hm)/sm
446        qm = qm + zpkdp(i, k)*(zq(i,k)-qm)/sm
447#ifdef ISO
[5158]448              DO ixt=1,ntraciso
[3927]449                 xtm(ixt)=xtm(ixt) &
450                    +zpkdp(i,k)*(zxt(ixt,i,k)-xtm(ixt))/sm
451              enddo
452#endif
453        down = .FALSE.
454        IF (k1/=limbas) THEN
455          IF (hm<zh(i,k1-1)) down = .TRUE.
456        END IF
457        IF (down) THEN
458          k1 = k1 - 1
459          k = k1
460        ELSE
461          IF ((k2==limhau)) GO TO 8021
462          IF ((zh(i,k2+1)>=hm)) GO TO 8021
463          k2 = k2 + 1
464          k = k2
465        END IF
466        GO TO 8020
4678021    CONTINUE
468        ! ------------ nouveau profil : constant (valeur moyenne)
469        DO k = k1, k2
470          zh(i, k) = hm
471          zq(i, k) = qm
472#ifdef ISO
[5158]473                DO ixt=1,ntraciso
[3927]474                   zxt(ixt,i,k)=xtm(ixt)
475                enddo
476#endif
477        END DO
478        k2 = k2 + 1
479      END IF
480      GO TO 8000
4818001  CONTINUE
482    END IF
483  END DO
484
485#ifdef ISO
486      ! cam verif
487#ifdef ISOVERIF
[5158]488      DO i=1,klon
489         DO k=limbas,limhau
490           DO ixt=1,ntraciso
[5103]491             CALL iso_verif_noNAN(zxt(ixt,i,k),'ajsec 428')
[3927]492           enddo !do ixt=1,niso           
[5117]493           IF (iso_eau.gt.0) THEN
[5103]494             CALL iso_verif_egalite_choix(zq(i,k),zxt(iso_eau,i,k), &
[3927]495                 'ajsec 432',errmax,errmaxrel)
[5116]496           endif !if (iso_eau.gt.0) THEN
[3927]497#ifdef ISOTRAC     
[5103]498           CALL iso_verif_traceur(zxt(1,i,k),'ajsec 436')
[3927]499#endif           
500         enddo !do k=limbas,limhau
501      enddo !do i=1,klon
502#endif
503      ! end cam verif
504#endif   
505
506
507  DO k = limbas, limhau
508    DO i = 1, klon
509      d_t(i, k) = zh(i, k)*zpk(i, k)/rcpd - t(i, k)
510      d_q(i, k) = zq(i, k) - q(i, k)
511#ifdef ISO
[5158]512         DO ixt=1,ntraciso
[3927]513            d_xt(ixt,i,k)=zxt(ixt,i,k)-xt(ixt,i,k)
514         enddo
515#endif
516    END DO
517  END DO
518
519  IF (limbas>1) THEN
520    DO k = 1, limbas - 1
521      DO i = 1, klon
522        d_t(i, k) = 0.0
523        d_q(i, k) = 0.0
524#ifdef ISO
[5158]525         DO ixt=1,ntraciso
[3927]526            d_xt(ixt,i,k)=0.0
527         enddo
528#endif
529      END DO
530    END DO
531  END IF
532
533  IF (limhau<klev) THEN
534    DO k = limhau + 1, klev
535      DO i = 1, klon
536        d_t(i, k) = 0.0
537        d_q(i, k) = 0.0
538#ifdef ISO
[5158]539         DO ixt=1,ntraciso
[3927]540            d_xt(ixt,i,k)=0.0
541         enddo
542#endif
543      END DO
544    END DO
545  END IF
546
547  IF (.NOT. mixq) THEN
548    DO k = 1, klev
549      DO i = 1, klon
550        d_q(i, k) = 0.0
551#ifdef ISO
[5158]552         DO ixt=1,ntraciso
[3927]553            d_xt(ixt,i,k)=0.0
554         enddo
555#endif
556      END DO
557    END DO
558  END IF
559
560
561#ifdef ISO
562      ! cam verif
563#ifdef ISOVERIF
[5158]564      DO i = 1, klon
565        DO k = limbas, limhau
[5117]566         IF (iso_eau.gt.0) THEN
[5103]567          CALL iso_verif_egalite(d_q(i,k),d_xt(iso_eau,i,k),'ajsec 270')
[3927]568         endif
569#ifdef ISOTRAC     
[5103]570        CALL iso_verif_traceur_justmass(d_xt(1,i,k),'physiq 3045')
[3927]571#endif         
572         enddo
573      enddo
574#endif
575      ! end cam verif
576#endif
577
[5105]578
[3927]579END SUBROUTINE ajsec_convv2
580SUBROUTINE ajsec_old(paprs, pplay, t, d_t)
581  USE dimphy
[5144]582        USE lmdz_yomcst
583
[3927]584  IMPLICIT 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  REAL paprs(klon, klev+1), pplay(klon, klev)
595  REAL t(klon, klev)
596  REAL d_t(klon, klev)
597
598  REAL local_h(klon, klev)
599  REAL hm, sm
600  LOGICAL modif(klon), down
601  INTEGER i, l, l1, l2
602  ! ------------------------------------- detection des profils a modifier
603  DO i = 1, klon
604    modif(i) = .FALSE.
605  END DO
606
607  DO l = 1, klev
608    DO i = 1, klon
609      local_h(i, l) = rcpd*t(i, l)/(pplay(i,l)**rkappa)
610    END DO
611  END DO
612
613  DO l = 2, klev
614    DO i = 1, klon
615      IF (local_h(i,l)<local_h(i,l-1)) THEN
616        modif(i) = .TRUE.
617      ELSE
618        modif(i) = modif(i)
619      END IF
620    END DO
621  END DO
622  ! ------------------------------------- correction des profils instables
623  DO i = 1, klon
624    IF (modif(i)) THEN
625      l2 = 1
6268000  CONTINUE
627      l2 = l2 + 1
628      IF (l2>klev) GO TO 8001
629      IF (local_h(i,l2)<local_h(i,l2-1)) THEN
630        l1 = l2 - 1
631        l = l1
632        sm = pplay(i, l2)**rkappa*(paprs(i,l2)-paprs(i,l2+1))
633        hm = local_h(i, l2)
6348020    CONTINUE
635        sm = sm + pplay(i, l)**rkappa*(paprs(i,l)-paprs(i,l+1))
636        hm = hm + pplay(i, l)**rkappa*(paprs(i,l)-paprs(i,l+1))*(local_h(i,l) &
637          -hm)/sm
638        down = .FALSE.
639        IF (l1/=1) THEN
640          IF (hm<local_h(i,l1-1)) THEN
641            down = .TRUE.
642          END IF
643        END IF
644        IF (down) THEN
645          l1 = l1 - 1
646          l = l1
647        ELSE
648          IF ((l2==klev)) GO TO 8021
649          IF ((local_h(i,l2+1)>=hm)) GO TO 8021
650          l2 = l2 + 1
651          l = l2
652        END IF
653        GO TO 8020
6548021    CONTINUE
655        ! ------------ nouveau profil : constant (valeur moyenne)
656        DO l = l1, l2
657          local_h(i, l) = hm
658        END DO
659        l2 = l2 + 1
660      END IF
661      GO TO 8000
6628001  CONTINUE
663    END IF
664  END DO
665
666  DO l = 1, klev
667    DO i = 1, klon
668      d_t(i, l) = local_h(i, l)*(pplay(i,l)**rkappa)/rcpd - t(i, l)
669    END DO
670  END DO
671
[5105]672
[3927]673END SUBROUTINE ajsec_old
Note: See TracBrowser for help on using the repository browser.