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

Last change on this file since 5271 was 5271, checked in by abarral, 26 hours ago

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