source: LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_thermcell_flux2.F90 @ 5209

Last change on this file since 5209 was 5158, checked in by abarral, 3 months ago

Add missing klon on strataer_emiss_mod.F90
Correct various missing explicit declarations
Replace tabs by spaces (tabs are not part of the fortran charset)
Continue cleaning modules
Removed unused arguments and variables
  • 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:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 17.2 KB
Line 
1MODULE lmdz_thermcell_flux2
2
3  ! $Id: lmdz_thermcell_flux2.F90 5158 2024-08-02 12:12:03Z fairhead $
4
5CONTAINS
6
7  SUBROUTINE thermcell_flux2(ngrid, nlay, ptimestep, masse, &
8          lalim, lmax, alim_star, &
9          entr_star, detr_star, f, rhobarz, zlev, zw2, fm, entr, &
10          detr, zqla, lev_out, lunout1, igout)
11    !IM 060508    &       detr,zqla,zmax,lev_out,lunout,igout)
12
13
14    !---------------------------------------------------------------------------
15    !thermcell_flux: deduction des flux
16    !---------------------------------------------------------------------------
17
18    USE lmdz_thermcell_ini, ONLY: prt_level, iflag_thermals_optflux
19    USE lmdz_abort_physic, ONLY: abort_physic
20    IMPLICIT NONE
21
22    ! arguments
23    INTEGER, INTENT(IN) :: ngrid, nlay
24    REAL, INTENT(IN) :: ptimestep
25    REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: masse
26    INTEGER, INTENT(IN), DIMENSION(ngrid) :: lalim, lmax
27    REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: alim_star, entr_star, detr_star
28    REAL, INTENT(IN), DIMENSION(ngrid) :: f
29    REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: rhobarz
30    REAL, INTENT(IN), DIMENSION(ngrid, nlay + 1) :: zw2, zlev
31    ! FH : laisser ca le temps de verifier qu'on a bien fait de commenter les
32    !      lignes faisant apparaitre zqla, zmax ...
33    !     REAL, INTENT(IN), DIMENSION(ngrid) :: zmax(ngrid)
34    !     enlever aussi zqla
35    REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: zqla  ! not used
36    INTEGER, INTENT(IN) :: lev_out, lunout1
37
38    REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: entr, detr
39    REAL, INTENT(OUT), DIMENSION(ngrid, nlay + 1) :: fm
40
41    ! local
42    INTEGER ig, l
43    INTEGER igout, lout
44    REAL zfm
45    INTEGER ncorecfm1, ncorecfm2, ncorecfm3, ncorecalpha
46    INTEGER ncorecfm4, ncorecfm5, ncorecfm6, ncorecfm7, ncorecfm8
47
48    REAL f_old, ddd0, eee0, ddd, eee, zzz
49
50    REAL, SAVE :: fomass_max = 0.5
51    REAL, SAVE :: alphamax = 0.7
52    !$OMP THREADPRIVATE(fomass_max,alphamax)
53
54    LOGICAL check_debug, labort_physic
55
56    CHARACTER (LEN = 20) :: modname = 'thermcell_flux2'
57    CHARACTER (LEN = 80) :: abort_message
58
59    ncorecfm1 = 0
60    ncorecfm2 = 0
61    ncorecfm3 = 0
62    ncorecfm4 = 0
63    ncorecfm5 = 0
64    ncorecfm6 = 0
65    ncorecfm7 = 0
66    ncorecfm8 = 0
67    ncorecalpha = 0
68
69    !initialisation
70    fm(:, :) = 0.
71
72    IF (prt_level>=10) THEN
73      WRITE(lunout1, *) 'Dans thermcell_flux 0'
74      WRITE(lunout1, *) 'flux base ', f(igout)
75      WRITE(lunout1, *) 'lmax ', lmax(igout)
76      WRITE(lunout1, *) 'lalim ', lalim(igout)
77      WRITE(lunout1, *) 'ig= ', igout
78      WRITE(lunout1, *) ' l E*    A*     D*  '
79      WRITE(lunout1, '(i4,3e15.5)') (l, entr_star(igout, l), alim_star(igout, l), detr_star(igout, l) &
80              , l = 1, lmax(igout))
81    endif
82
83
84    !-------------------------------------------------------------------------
85    ! Verification de la nullite des entrainement et detrainement au dessus
86    ! de lmax(ig)
87    ! Active uniquement si check_debug=.TRUE. ou prt_level>=10
88    !-------------------------------------------------------------------------
89
90    check_debug = .false..or.prt_level>=10
91
92    IF (check_debug) THEN
93      DO l = 1, nlay
94        DO ig = 1, ngrid
95          IF (l<=lmax(ig)) THEN
96            IF (entr_star(ig, l)>1.) THEN
97              PRINT*, 'WARNING thermcell_flux 1 ig,l,lmax(ig)', ig, l, lmax(ig)
98              PRINT*, 'entr_star(ig,l)', entr_star(ig, l)
99              PRINT*, 'alim_star(ig,l)', alim_star(ig, l)
100              PRINT*, 'detr_star(ig,l)', detr_star(ig, l)
101            endif
102          else
103            IF (abs(entr_star(ig, l)) + abs(alim_star(ig, l)) + abs(detr_star(ig, l))>0.) THEN
104              PRINT*, 'cas 1 : ig,l,lmax(ig)', ig, l, lmax(ig)
105              PRINT*, 'entr_star(ig,l)', entr_star(ig, l)
106              PRINT*, 'alim_star(ig,l)', alim_star(ig, l)
107              PRINT*, 'detr_star(ig,l)', detr_star(ig, l)
108              abort_message = ''
109              labort_physic = .TRUE.
110              CALL abort_physic (modname, abort_message, 1)
111            endif
112          endif
113        enddo
114      enddo
115    endif
116
117    !-------------------------------------------------------------------------
118    ! Multiplication par le flux de masse issu de la femreture
119    !-------------------------------------------------------------------------
120
121    DO l = 1, nlay
122      entr(:, l) = f(:) * (entr_star(:, l) + alim_star(:, l))
123      detr(:, l) = f(:) * detr_star(:, l)
124    enddo
125
126    IF (prt_level>=10) THEN
127      WRITE(lunout1, *) 'Dans thermcell_flux 1'
128      WRITE(lunout1, *) 'flux base ', f(igout)
129      WRITE(lunout1, *) 'lmax ', lmax(igout)
130      WRITE(lunout1, *) 'lalim ', lalim(igout)
131      WRITE(lunout1, *) 'ig= ', igout
132      WRITE(lunout1, *) ' l   E    D     W2'
133      WRITE(lunout1, '(i4,3e15.5)') (l, entr(igout, l), detr(igout, l) &
134              , zw2(igout, l + 1), l = 1, lmax(igout))
135    endif
136
137    fm(:, 1) = 0.
138    DO l = 1, nlay
139      DO ig = 1, ngrid
140        IF (l<lmax(ig)) THEN
141          fm(ig, l + 1) = fm(ig, l) + entr(ig, l) - detr(ig, l)
142        elseif(l==lmax(ig)) THEN
143          fm(ig, l + 1) = 0.
144          detr(ig, l) = fm(ig, l) + entr(ig, l)
145        else
146          fm(ig, l + 1) = 0.
147        endif
148      enddo
149    enddo
150
151
152
153    ! Test provisoire : pour comprendre pourquoi on corrige plein de fois
154    ! le cas fm6, on commence par regarder une premiere fois avant les
155    ! autres corrections.
156
157    DO l = 1, nlay
158      DO ig = 1, ngrid
159        IF (detr(ig, l)>fm(ig, l)) THEN
160          ncorecfm8 = ncorecfm8 + 1
161          !              igout=ig
162        endif
163      enddo
164    enddo
165
166    !      if (prt_level.ge.10) &
167    !    &    CALL printflux(ngrid,nlay,lunout1,igout,f,lmax,lalim, &
168    !    &    ptimestep,masse,entr,detr,fm,'2  ')
169
170
171
172    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
173    ! FH Version en cours de test;
174    ! par rapport a thermcell_flux, on fait une grande boucle sur "l"
175    ! et on modifie le flux avec tous les contrôles appliques d'affilee
176    ! pour la meme couche
177    ! Momentanement, on duplique le calcule du flux pour pouvoir comparer
178    ! les flux avant et apres modif
179    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
180
181    DO l = 1, nlay
182
183      DO ig = 1, ngrid
184        IF (l<lmax(ig)) THEN
185          fm(ig, l + 1) = fm(ig, l) + entr(ig, l) - detr(ig, l)
186        elseif(l==lmax(ig)) THEN
187          fm(ig, l + 1) = 0.
188          detr(ig, l) = fm(ig, l) + entr(ig, l)
189        else
190          fm(ig, l + 1) = 0.
191        endif
192      enddo
193
194
195      !-------------------------------------------------------------------------
196      ! Verification de la positivite des flux de masse
197      !-------------------------------------------------------------------------
198
199      !     do l=1,nlay
200      DO ig = 1, ngrid
201        IF (fm(ig, l + 1)<0.) THEN
202          !              PRINT*,'fm1<0',l+1,lmax(ig),fm(ig,l+1)
203          ncorecfm1 = ncorecfm1 + 1
204          fm(ig, l + 1) = fm(ig, l)
205          detr(ig, l) = entr(ig, l)
206        endif
207      enddo
208      !     enddo
209
210      IF (prt_level>=10) &
211              WRITE(lunout1, '(i4,4e14.4)') l, masse(igout, l) / ptimestep, &
212                      entr(igout, l), detr(igout, l), fm(igout, l + 1)
213
214      !-------------------------------------------------------------------------
215      !Test sur fraca croissant
216      !-------------------------------------------------------------------------
217      IF (iflag_thermals_optflux==0) THEN
218        !     do l=1,nlay
219        DO ig = 1, ngrid
220          IF (l>=lalim(ig).AND.l<=lmax(ig) &
221                  .AND.(zw2(ig, l + 1)>1.e-10).AND.(zw2(ig, l)>1.e-10)) THEN
222            !  zzz est le flux en l+1 a frac constant
223            zzz = fm(ig, l) * rhobarz(ig, l + 1) * zw2(ig, l + 1)  &
224                    / (rhobarz(ig, l) * zw2(ig, l))
225            IF (fm(ig, l + 1)>zzz) THEN
226              detr(ig, l) = detr(ig, l) + fm(ig, l + 1) - zzz
227              fm(ig, l + 1) = zzz
228              ncorecfm4 = ncorecfm4 + 1
229            endif
230          endif
231        enddo
232        !     enddo
233      endif
234
235      IF (prt_level>=10) &
236              WRITE(lunout1, '(i4,4e14.4)') l, masse(igout, l) / ptimestep, &
237                      entr(igout, l), detr(igout, l), fm(igout, l + 1)
238
239
240      !-------------------------------------------------------------------------
241      !test sur flux de masse croissant
242      !-------------------------------------------------------------------------
243      IF (iflag_thermals_optflux==0) THEN
244        !     do l=1,nlay
245        DO ig = 1, ngrid
246          IF ((fm(ig, l + 1)>fm(ig, l)).AND.(l>lalim(ig))) THEN
247            f_old = fm(ig, l + 1)
248            fm(ig, l + 1) = fm(ig, l)
249            detr(ig, l) = detr(ig, l) + f_old - fm(ig, l + 1)
250            ncorecfm5 = ncorecfm5 + 1
251          endif
252        enddo
253        !     enddo
254      endif
255
256      IF (prt_level>=10) &
257              WRITE(lunout1, '(i4,4e14.4)') l, masse(igout, l) / ptimestep, &
258                      entr(igout, l), detr(igout, l), fm(igout, l + 1)
259
260      !fin 1.EQ.0
261      !-------------------------------------------------------------------------
262      !detr ne peut pas etre superieur a fm
263      !-------------------------------------------------------------------------
264
265      IF(1==1) THEN
266        !     do l=1,nlay
267
268        labort_physic = .FALSE.
269        DO ig = 1, ngrid
270          IF (entr(ig, l)<0.) THEN
271            labort_physic = .TRUE.
272            igout = ig
273            lout = l
274          endif
275        enddo
276
277        IF (labort_physic) THEN
278          PRINT*, 'N1 ig,l,entr', igout, lout, entr(igout, lout)
279          abort_message = 'entr negatif'
280          CALL abort_physic (modname, abort_message, 1)
281        endif
282
283        DO ig = 1, ngrid
284          IF (detr(ig, l)>fm(ig, l)) THEN
285            ncorecfm6 = ncorecfm6 + 1
286            detr(ig, l) = fm(ig, l)
287            entr(ig, l) = fm(ig, l + 1)
288
289            ! Dans le cas ou on est au dessus de la couche d'alimentation et que le
290            ! detrainement est plus fort que le flux de masse, on stope le thermique.
291            !test:on commente
292            !               if (l.gt.lalim(ig)) THEN
293            !                  lmax(ig)=l
294            !                  fm(ig,l+1)=0.
295            !                  entr(ig,l)=0.
296            !               else
297            !                  ncorecfm7=ncorecfm7+1
298            !               endif
299          endif
300
301          IF(l>lmax(ig)) THEN
302            detr(ig, l) = 0.
303            fm(ig, l + 1) = 0.
304            entr(ig, l) = 0.
305          endif
306        enddo
307
308        labort_physic = .FALSE.
309        DO ig = 1, ngrid
310          IF (entr(ig, l)<0.) THEN
311            labort_physic = .TRUE.
312            igout = ig
313          endif
314        enddo
315        IF (labort_physic) THEN
316          ig = igout
317          PRINT*, 'ig,l,lmax(ig)', ig, l, lmax(ig)
318          PRINT*, 'entr(ig,l)', entr(ig, l)
319          PRINT*, 'fm(ig,l)', fm(ig, l)
320          abort_message = 'probleme dans thermcell flux'
321          CALL abort_physic (modname, abort_message, 1)
322        endif
323
324
325        !     enddo
326      endif
327
328      IF (prt_level>=10) &
329              WRITE(lunout1, '(i4,4e14.4)') l, masse(igout, l) / ptimestep, &
330                      entr(igout, l), detr(igout, l), fm(igout, l + 1)
331
332      !-------------------------------------------------------------------------
333      !fm ne peut pas etre negatif
334      !-------------------------------------------------------------------------
335
336      !     do l=1,nlay
337      DO ig = 1, ngrid
338        IF (fm(ig, l + 1)<0.) THEN
339          detr(ig, l) = detr(ig, l) + fm(ig, l + 1)
340          fm(ig, l + 1) = 0.
341          ncorecfm2 = ncorecfm2 + 1
342        endif
343      enddo
344
345      labort_physic = .FALSE.
346      DO ig = 1, ngrid
347        IF (detr(ig, l)<0.) THEN
348          labort_physic = .TRUE.
349          igout = ig
350        endif
351      enddo
352      IF (labort_physic) THEN
353        ig = igout
354        PRINT*, 'cas 2 : ig,l,lmax(ig)', ig, l, lmax(ig)
355        PRINT*, 'detr(ig,l)', detr(ig, l)
356        PRINT*, 'fm(ig,l)', fm(ig, l)
357        abort_message = 'probleme dans thermcell flux'
358        CALL abort_physic (modname, abort_message, 1)
359      endif
360      !    enddo
361
362      IF (prt_level>=10) &
363              WRITE(lunout1, '(i4,4e14.4)') l, masse(igout, l) / ptimestep, &
364                      entr(igout, l), detr(igout, l), fm(igout, l + 1)
365
366      !-----------------------------------------------------------------------
367      !la fraction couverte ne peut pas etre superieure a 1
368      !-----------------------------------------------------------------------
369
370
371      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
372      ! FH Partie a revisiter.
373      ! Il semble qu'etaient codees ici deux optiques dans le cas
374      ! F/ (rho *w) > 1
375      ! soit limiter la hauteur du thermique en considerant que c'est
376      ! la derniere chouche, soit limiter F a rho w.
377      ! Dans le second cas, il faut en fait limiter a un peu moins
378      ! que ca parce qu'on a des 1 / ( 1 -alpha) un peu plus loin
379      ! dans thermcell_main et qu'il semble de toutes facons deraisonable
380      ! d'avoir des fractions de 1..
381      ! Ci dessous, et dans l'etat actuel, le premier des  deux if est
382      ! sans doute inutile.
383      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
384
385      !    do l=1,nlay
386      DO ig = 1, ngrid
387        IF (zw2(ig, l + 1)>1.e-10) THEN
388          zfm = rhobarz(ig, l + 1) * zw2(ig, l + 1) * alphamax
389          IF (fm(ig, l + 1) > zfm) THEN
390            f_old = fm(ig, l + 1)
391            fm(ig, l + 1) = zfm
392            detr(ig, l) = detr(ig, l) + f_old - fm(ig, l + 1)
393            !             lmax(ig)=l+1
394            !             zmax(ig)=zlev(ig,lmax(ig))
395            !             PRINT*,'alpha>1',l+1,lmax(ig)
396            ncorecalpha = ncorecalpha + 1
397          endif
398        endif
399      enddo
400      !    enddo
401
402      IF (prt_level>=10) &
403              WRITE(lunout1, '(i4,4e14.4)') l, masse(igout, l) / ptimestep, &
404                      entr(igout, l), detr(igout, l), fm(igout, l + 1)
405
406      ! Fin de la grande boucle sur les niveaux verticaux
407    enddo
408
409    !      if (prt_level.ge.10) &
410    !    &    CALL printflux(ngrid,nlay,lunout1,igout,f,lmax,lalim, &
411    !    &    ptimestep,masse,entr,detr,fm,'8  ')
412
413
414    !-----------------------------------------------------------------------
415    ! On fait en sorte que la quantite totale d'air entraine dans le
416    ! panache ne soit pas trop grande comparee a la masse de la maille
417    !-----------------------------------------------------------------------
418
419    IF (1==1) THEN
420      labort_physic = .FALSE.
421      DO l = 1, nlay - 1
422        DO ig = 1, ngrid
423          eee0 = entr(ig, l)
424          ddd0 = detr(ig, l)
425          eee = entr(ig, l) - masse(ig, l) * fomass_max / ptimestep
426          ddd = detr(ig, l) - eee
427          IF (eee>0.) THEN
428            ncorecfm3 = ncorecfm3 + 1
429            entr(ig, l) = entr(ig, l) - eee
430            IF (ddd>0.) THEN
431              !   l'entrainement est trop fort mais l'exces peut etre compense par une
432              !   diminution du detrainement)
433              detr(ig, l) = ddd
434            else
435              !   l'entrainement est trop fort mais l'exces doit etre compense en partie
436              !   par un entrainement plus fort dans la couche superieure
437              IF(l==lmax(ig)) THEN
438                detr(ig, l) = fm(ig, l) + entr(ig, l)
439              else
440                IF(l>=lmax(ig).AND.0==1) THEN
441                  igout = ig
442                  lout = l
443                  labort_physic = .TRUE.
444                endif
445                entr(ig, l + 1) = entr(ig, l + 1) - ddd
446                detr(ig, l) = 0.
447                fm(ig, l + 1) = fm(ig, l) + entr(ig, l)
448                detr(ig, l) = 0.
449              endif
450            endif
451          endif
452        enddo
453      enddo
454      IF (labort_physic) THEN
455        ig = igout
456        l = lout
457        PRINT*, 'ig,l', ig, l
458        PRINT*, 'eee0', eee0
459        PRINT*, 'ddd0', ddd0
460        PRINT*, 'eee', eee
461        PRINT*, 'ddd', ddd
462        PRINT*, 'entr', entr(ig, l)
463        PRINT*, 'detr', detr(ig, l)
464        PRINT*, 'masse', masse(ig, l)
465        PRINT*, 'fomass_max', fomass_max
466        PRINT*, 'masse(ig,l)*fomass_max/ptimestep', masse(ig, l) * fomass_max / ptimestep
467        PRINT*, 'ptimestep', ptimestep
468        PRINT*, 'lmax(ig)', lmax(ig)
469        PRINT*, 'fm(ig,l+1)', fm(ig, l + 1)
470        PRINT*, 'fm(ig,l)', fm(ig, l)
471        abort_message = 'probleme dans thermcell_flux'
472        CALL abort_physic (modname, abort_message, 1)
473      endif
474    endif
475
476    !              ddd=detr(ig)-entre
477    !on s assure que tout s annule bien en zmax
478    DO ig = 1, ngrid
479      fm(ig, lmax(ig) + 1) = 0.
480      entr(ig, lmax(ig)) = 0.
481      detr(ig, lmax(ig)) = fm(ig, lmax(ig)) + entr(ig, lmax(ig))
482    enddo
483
484    !-----------------------------------------------------------------------
485    ! Impression du nombre de bidouilles qui ont ete necessaires
486    !-----------------------------------------------------------------------
487
488    !IM 090508 beg
489    !     if (ncorecfm1+ncorecfm2+ncorecfm3+ncorecfm4+ncorecfm5+ncorecalpha > 0 ) THEN
490    !         PRINT*,'PB thermcell : on a du coriger ',ncorecfm1,'x fm1',&
491    !   &     ncorecfm2,'x fm2',ncorecfm3,'x fm3 et', &
492    !   &     ncorecfm4,'x fm4',ncorecfm5,'x fm5 et', &
493    !   &     ncorecfm6,'x fm6', &
494    !   &     ncorecfm7,'x fm7', &
495    !   &     ncorecfm8,'x fm8', &
496    !   &     ncorecalpha,'x alpha'
497    !     endif
498    !IM 090508 end
499
500    !      if (prt_level.ge.10) &
501    !    &    CALL printflux(ngrid,nlay,lunout1,igout,f,lmax,lalim, &
502    !    &    ptimestep,masse,entr,detr,fm,'fin')
503
504    RETURN
505  END
506END MODULE lmdz_thermcell_flux2
Note: See TracBrowser for help on using the repository browser.