source: LMDZ6/branches/LMDZ-tracers/libf/dyn3dmem/qminimum_loc.F @ 4007

Last change on this file since 4007 was 3891, checked in by dcugnet, 4 years ago
  • Bugs corrections:
    • sequential gcm fixed
    • parallel gcm compilation fixed ; to be tested
  • Some generic operations moved from infotrac to readTracFile
  • Fixed algebrical reduction routine, used in the isotopes parameters file.
  • Additional component "comp" in the tracers descriptor derived type "tra",

specifying the model component name(s) (cf. tracers sections) it belongs.

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