source: LMDZ6/branches/Amaury_dev/libf/dyn3dmem/qminimum_loc.f90 @ 5123

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

Correct various minor mistakes from previous commits
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:keywords set to Id
File size: 10.0 KB
Line 
1! $Id: qminimum_loc.f90 5123 2024-07-25 06:45:50Z abarral $
2
3SUBROUTINE qminimum_loc(q, nqtot, deltap)
4  USE parallel_lmdz
5  USE infotrac, ONLY: niso, ntiso, iqIsoPha, tracers, &
6          isoCheck, min_qParent
7  USE lmdz_strings, ONLY: strIdx
8  USE lmdz_readTracFiles, ONLY: addPhase
9  USE lmdz_iniprint, ONLY: lunout, prt_level
10  IMPLICIT none
11  !
12  !  -- Objet : Traiter les valeurs trop petites (meme negatives)
13  !         pour l'eau vapeur et l'eau liquide
14  !
15  include "dimensions.h"
16  include "paramet.h"
17  !
18  INTEGER :: nqtot ! CRisi: on remplace nq par nqtot
19  REAL :: q(ijb_u:ije_u, llm, nqtot), deltap(ijb_u:ije_u, llm)
20  !
21  LOGICAL, SAVE :: first = .TRUE.
22  INTEGER, SAVE :: iq_vap, iq_liq        ! indices pour l'eau vapeur/liquide
23  !$OMP THREADPRIVATE(iq_vap, iq_liq, first)
24  REAL, PARAMETER :: seuil_vap = 1.0e-10 ! seuil pour l'eau vapeur
25  REAL, PARAMETER :: seuil_liq = 1.0e-11 ! seuil pour l'eau liquide
26  !
27  !  NB. ....( Il est souhaitable mais non obligatoire que les valeurs des
28  !        parametres seuil_vap, seuil_liq soient pareilles a celles
29  !        qui  sont utilisees dans la routine    ADDFI       )
30  ! .................................................................
31  !
32  !DC iq_val and iq_liq are usable for q only, NOT for q_follow
33  !   and zx_defau_diag (crash if iq_val/liq==3) => vapor/liquid
34  !   water at hardcoded indices 1/2 in these variables
35  INTEGER :: i, k, iq
36  REAL :: zx_defau, zx_abc, zx_pump(ijb_u:ije_u), pompe
37
38  REAL :: zx_defau_diag(ijb_u:ije_u, llm, 2)
39  REAL :: q_follow(ijb_u:ije_u, llm, 2)
40  !
41  INTEGER :: imprim
42  SAVE imprim
43  DATA imprim /0/
44  !$OMP THREADPRIVATE(imprim)
45  INTEGER :: ijb, ije
46  INTEGER :: Index_pump(ij_end - ij_begin + 1)
47  INTEGER :: nb_pump
48  INTEGER :: ixt
49  INTEGER :: iso_verif_noNaN_nostop
50
51  !$OMP BARRIER
52
53  !WRITE(lunout,*) 'qminimum 52: entree'
54  IF(first) THEN
55    iq_vap = strIdx(tracers(:)%name, addPhase('H2O', 'g'))
56    iq_liq = strIdx(tracers(:)%name, addPhase('H2O', 'l'))
57    first = .FALSE.
58  END IF
59  !
60  ! Quand l'eau liquide est trop petite (ou negative), on prend
61  ! l'eau vapeur de la meme couche et la convertit en eau liquide
62  ! (sans changer la temperature !)
63  !
64
65  CALL check_isotopes(q, ij_begin, ij_end, 'qminimum 52')
66
67  ijb = ij_begin
68  ije = ij_end
69
70  DO k = 1, llm
71    !$OMP DO SCHEDULE(STATIC)
72    DO i = ijb, ije
73      zx_defau_diag(i, k, 1) = 0.0
74      zx_defau_diag(i, k, 2) = 0.0
75      q_follow(i, k, 1) = q(i, k, iq_vap)
76      q_follow(i, k, 2) = q(i, k, iq_liq)
77    ENDDO
78    !$OMP END DO NOWAIT
79  ENDDO
80
81  !WRITE(lunout,*) 'qminimum 57'
82  DO k = 1, llm
83    !$OMP DO SCHEDULE(STATIC)
84    DO i = ijb, ije
85      IF (seuil_liq - q(i, k, iq_liq) > 0.d0) THEN
86        IF (niso > 0) zx_defau_diag(i, k, 2) = AMAX1 &
87                (seuil_liq - q(i, k, iq_liq), 0.0)
88
89        q(i, k, iq_vap) = q(i, k, iq_vap) + q(i, k, iq_liq) - seuil_liq
90        q(i, k, iq_liq) = seuil_liq
91      endif
92    END DO
93    !$OMP END DO NOWAIT
94  END DO
95
96  !
97  ! Quand l'eau vapeur est trop faible (ou negative), on complete
98  ! le defaut en prennant de l'eau vapeur de la couche au-dessous.
99  !
100  !WRITE(lunout,*) 'qminimum 81'
101  DO k = llm, 2, -1
102    !cc      zx_abc = dpres(k) / dpres(k-1)
103    !$OMP DO SCHEDULE(STATIC)
104    DO i = ijb, ije
105
106      IF (seuil_vap - q(i, k, iq_vap) > 0.d0) THEN
107        IF (niso > 0) zx_defau_diag(i, k, 1) &
108                = AMAX1(seuil_vap - q(i, k, iq_vap), 0.0)
109
110        q(i, k - 1, iq_vap) = q(i, k - 1, iq_vap) - (seuil_vap &
111                - q(i, k, iq_vap)) * deltap(i, k) / deltap(i, k - 1)
112        q(i, k, iq_vap) = seuil_vap
113
114      endif
115    ENDDO
116    !$OMP END DO NOWAIT
117  ENDDO
118
119  !
120  ! Quand il s'agit de la premiere couche au-dessus du sol, on
121  ! doit imprimer un message d'avertissement (saturation possible).
122  !
123  !WRITE(lunout,*) 'qminimum 106'
124  nb_pump = 0
125  !$OMP DO SCHEDULE(STATIC)
126  DO i = ijb, ije
127    zx_pump(i) = AMAX1(0.0, seuil_vap - q(i, 1, iq_vap))
128    q(i, 1, iq_vap) = AMAX1(q(i, 1, iq_vap), seuil_vap)
129    IF (zx_pump(i) > 0.0) THEN
130      nb_pump = nb_pump + 1
131      Index_pump(nb_pump) = i
132    ENDIF
133  ENDDO
134  !$OMP END DO NOWAIT
135  ! pompe = SSUM(ije-ijb+1,zx_pump(ijb),1)
136
137  IF (imprim<=100 .AND. nb_pump > 0) THEN
138    PRINT *, 'ATT!:on pompe de l eau au sol'
139    DO i = 1, nb_pump
140      imprim = imprim + 1
141      PRINT*, '  en ', index_pump(i), zx_pump(index_pump(i))
142    ENDDO
143  ENDIF
144
145  !WRITE(lunout,*) 'qminimum 128'
146  IF (niso > 0) THEN
147    !WRITE(lunout,*) 'qminimum 140'
148    ! CRisi: traiter de même les traceurs d'eau
149    ! Mais il faut les prendre à l'envers pour essayer de conserver la
150    ! masse.
151    ! 1) pompage dans le sol
152    ! On suppose que ce pompage se fait sans isotopes -> on ne modifie
153    ! rien ici et on croise les doigts pour que ça ne soit pas trop
154    ! génant
155    ! en fait, si, c'est genant quand les isotopes doivent eux même transporter des
156    ! traceurs -> apporter aussi un peu d'isotopes... Combien?
157    ! Essayer tnat/2 = -500 permil? C'est déjà mieux que -1000
158    ! permil...
159    ! pb: que faire pour les traceurs?
160    !$OMP DO SCHEDULE(STATIC)
161    DO i = ijb, ije
162      IF (zx_pump(i)>0.0) THEN
163        q_follow(i, 1, 1) = q_follow(i, 1, 1) + zx_pump(i)
164      endif !if (zx_pump(i).gt.0.0) THEN
165    enddo !DO i = ijb, ije
166    !$OMP END DO NOWAIT
167
168    ! 2) transfert de vap vers les couches plus hautes
169    !WRITE(lunout,*) 'qminimum 158'
170    do k = 2, llm
171      !$OMP DO SCHEDULE(STATIC)
172      DO i = ijb, ije
173        IF (zx_defau_diag(i, k, 1)>0.0) THEN
174          ! on ajoute la vapeur en k
175          !  WRITE(lunout,*) 'i,k,q_follow(i,k-1,ivap)=',
176          ! :                 i,k,q_follow(i,k-1,1)
177          IF (q_follow(i, k - 1, 1)<min_qParent) THEN
178            WRITE(lunout, *) 'tmp qmin: on stoppe'
179            WRITE(lunout, *) 'zx_pump(i)=', zx_pump(i)
180            WRITE(lunout, *) 'q_follow(i,:,ivap)=', &
181                    q_follow(i, :, 1)
182            WRITE(lunout, *) 'k=', k
183            CALL abort_gcm("qminimum", "not enough vapor", 1)
184          endif
185          do ixt = 1, ntiso
186            ! WRITE(lunout,*) 'qmin 168: ixt=',ixt
187            ! WRITE(lunout,*) 'q(i,k,iqIsoPha(ixt,iq_vap))=',
188            ! :             q(i,k,iqIsoPha(ixt,iq_vap))
189            !            WRITE(lunout,*) 'zx_defau_diag(i,k,ivap)=',
190            ! :                  zx_defau_diag(i,k,1)
191            !            WRITE(lunout,*) 'q(i,k-1,iqIsoPha(ixt,iq_vap))=',
192            ! :                   q(i,k-1,iqIsoPha(ixt,iq_vap))
193
194            q(i, k, iqIsoPha(ixt, iq_vap)) = q(i, k, iqIsoPha(ixt, iq_vap)) &
195                    + zx_defau_diag(i, k, 1) &
196                            * q(i, k - 1, iqIsoPha(ixt, iq_vap)) / q_follow(i, k - 1, 1)
197
198            IF (isoCheck) THEN
199              IF(iso_verif_noNaN_nostop(q(i, k, iqIsoPha(ixt, iq_vap)), &
200                      'qminimum 155')==1) THEN
201                WRITE(*, *) 'i,k,ixt=', i, k, ixt
202                WRITE(*, *) 'q_follow(i,k-1,ivap)=', &
203                        q_follow(i, k - 1, 1)
204                WRITE(*, *) 'q(i,k,iqIsoPha(ixt,iq_vap))=', &
205                        q(i, k, iqIsoPha(ixt, iq_vap))
206                WRITE(*, *) 'zx_defau_diag(i,k,ivap)=', &
207                        zx_defau_diag(i, k, 1)
208                WRITE(*, *) 'q(i,k-1,iqIsoPha(ixt,iq_vap))=', &
209                        q(i, k - 1, iqIsoPha(ixt, iq_vap))
210                CALL abort_gcm("qminimum_loc", "stopped", 1)
211              endif
212            endif
213
214            ! et on la retranche en k-1
215            q(i, k - 1, iqIsoPha(ixt, iq_vap)) = &
216                    q(i, k - 1, iqIsoPha(ixt, iq_vap)) &
217                            - zx_defau_diag(i, k, 1) &
218                            * deltap(i, k) / deltap(i, k - 1) &
219                            * q(i, k - 1, iqIsoPha(ixt, iq_vap)) &
220                            / q_follow(i, k - 1, 1)
221
222            IF (isoCheck) THEN
223              IF (iso_verif_noNaN_nostop(&
224                      q(i, k - 1, iqIsoPha(ixt, iq_vap)), &
225                      'qminimum 175')==1) THEN
226                WRITE(*, *) 'k,i,ixt=', k, i, ixt
227                WRITE(*, *) 'q_follow(i,k-1,ivap)=', &
228                        q_follow(i, k - 1, 1)
229                WRITE(*, *) 'q(i,k,iqIsoPha(ixt,iq_vap))=', &
230                        q(i, k, iqIsoPha(ixt, iq_vap))
231                WRITE(*, *) 'zx_defau_diag(i,k,ivap)=', &
232                        zx_defau_diag(i, k, 1)
233                WRITE(*, *) 'q(i,k-1,iqIsoPha(ixt,iq_vap))=', &
234                        q(i, k - 1, iqIsoPha(ixt, iq_vap))
235                CALL abort_gcm("qminimum_loc", "stopped", 1)
236              endif
237            endif
238
239          enddo !do ixt=1,niso
240          q_follow(i, k, 1) = q_follow(i, k, 1) &
241                  + zx_defau_diag(i, k, 1)
242          q_follow(i, k - 1, 1) = q_follow(i, k - 1, 1) &
243                  - zx_defau_diag(i, k, 1) &
244                          * deltap(i, k) / deltap(i, k - 1)
245        endif !if (zx_defau_diag(i,k,1).gt.0.0) THEN
246      enddo !DO i = 1, ip1jmp1
247      !$OMP END DO NOWAIT
248    enddo !do k=2,llm
249
250    CALL check_isotopes(q, ijb, ije, 'qminimum 168')
251
252
253    ! 3) transfert d'eau de la vapeur au liquide
254    !WRITE(*,*) 'qminimum 164'
255    do k = 1, llm
256      !$OMP DO SCHEDULE(STATIC)
257      DO i = ijb, ije
258        IF (zx_defau_diag(i, k, 2)>0.0) THEN
259          ! on ajoute eau liquide en k en k
260          do ixt = 1, ntiso
261            q(i, k, iqIsoPha(ixt, iq_liq)) = q(i, k, iqIsoPha(ixt, iq_liq)) &
262                    + zx_defau_diag(i, k, 2) &
263                            * q(i, k, iqIsoPha(ixt, iq_vap)) / q_follow(i, k, 1)
264            ! et on la retranche à la vapeur en k
265            q(i, k, iqIsoPha(ixt, iq_vap)) = q(i, k, iqIsoPha(ixt, iq_vap)) &
266                    - zx_defau_diag(i, k, 2) &
267                            * q(i, k, iqIsoPha(ixt, iq_vap)) / q_follow(i, k, 1)
268          enddo !do ixt=1,niso
269          q_follow(i, k, 2) = q_follow(i, k, 2) &
270                  + zx_defau_diag(i, k, 2)
271          q_follow(i, k, 1) = q_follow(i, k, 1) &
272                  - zx_defau_diag(i, k, 2)
273        endif !if (zx_defau_diag(i,k,1).gt.0.0) THEN
274      enddo !DO i = ijb, ije
275      !$OMP END DO NOWAIT
276    enddo !do k=2,llm
277
278    CALL check_isotopes(q, ijb, ije, 'qminimum 197')
279
280  ENDIF !if (niso > 0) THEN
281  !WRITE(*,*) 'qminimum 188'
282  !$OMP BARRIER
283
284  !
285
286END SUBROUTINE qminimum_loc
Note: See TracBrowser for help on using the repository browser.