source: LMDZ6/branches/Ocean_skin/libf/dyn3dmem/qminimum_loc.F

Last change on this file was 4368, checked in by lguez, 19 months ago

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