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

Last change on this file since 5159 was 5159, checked in by abarral, 7 weeks ago

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