source: LMDZ6/branches/cirrus/libf/dyn3dmem/qminimum_loc.F @ 5210

Last change on this file since 5210 was 5203, checked in by Laurent Fairhead, 30 hours ago

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