source: LMDZ6/trunk/libf/dyn3dmem/qminimum_loc.F @ 3800

Last change on this file since 3800 was 3800, checked in by Laurent Fairhead, 4 years ago

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