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

Last change on this file since 5259 was 5246, checked in by abarral, 15 months ago

Convert fixed-form to free-form sources .F -> .{f,F}90
(WIP: some .F remain, will be handled in subsequent commits)

  • 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 5246 2024-10-21 12:58:45Z 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  IMPLICIT none
10  !
11  !  -- Objet : Traiter les valeurs trop petites (meme negatives)
12  !         pour l'eau vapeur et l'eau liquide
13  !
14  include "dimensions.h"
15  include "paramet.h"
16  include "iniprint.h"
17  !
18  INTEGER :: nqtot ! CRisi: on remplace nq par nqtot
19  REAL :: q(ijb_u:ije_u,llm,nqtot), deltap(ijb_u:ije_u,llm)
20  !
21  LOGICAL, SAVE :: first=.TRUE.
22  INTEGER, SAVE :: iq_vap, iq_liq        ! indices pour l'eau vapeur/liquide
23!$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
26  !
27  !  NB. ....( Il est souhaitable mais non obligatoire que les valeurs des
28  !        parametres seuil_vap, seuil_liq soient pareilles a celles
29  !        qui  sont utilisees dans la routine    ADDFI       )
30  ! .................................................................
31  !
32  !DC iq_val and iq_liq are usable for q only, NOT for q_follow
33  !   and zx_defau_diag (crash if iq_val/liq==3) => vapor/liquid
34  !   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)
40  !
41  REAL :: SSUM
42  EXTERNAL SSUM
43  !
44  INTEGER :: imprim
45  SAVE imprim
46  DATA imprim /0/
47!$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
54!$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
62  !
63  ! Quand l'eau liquide est trop petite (ou negative), on prend
64  ! l'eau vapeur de la meme couche et la convertit en eau liquide
65  ! (sans changer la temperature !)
66  !
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
74!$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
81!$OMP END DO NOWAIT
82  ENDDO
83
84  ! !write(lunout,*) 'qminimum 57'
85  DO k = 1, llm
86!$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
97!$OMP END DO NOWAIT
98  END DO
99
100  !
101  ! Quand l'eau vapeur est trop faible (ou negative), on complete
102  ! le defaut en prennant de l'eau vapeur de la couche au-dessous.
103  !
104  ! !write(lunout,*) 'qminimum 81'
105  DO k = llm, 2, -1
106  !cc      zx_abc = dpres(k) / dpres(k-1)
107!$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
121!$OMP END DO NOWAIT
122  ENDDO
123
124  !
125  ! Quand il s'agit de la premiere couche au-dessus du sol, on
126  ! doit imprimer un message d'avertissement (saturation possible).
127  !
128  ! !write(lunout,*) 'qminimum 106'
129  nb_pump=0
130!$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
139!$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?
165!$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
171!$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
176!$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
252!$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
261!$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
281!$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'
288!$OMP BARRIER
289
290  !
291  RETURN
292END SUBROUTINE qminimum_loc
Note: See TracBrowser for help on using the repository browser.