source: LMDZ6/trunk/libf/dyn3dmem/qminimum_loc.f90 @ 5279

Last change on this file since 5279 was 5272, checked in by abarral, 33 hours ago

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