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

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

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