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

Last change on this file since 5135 was 5119, checked in by abarral, 16 months ago

enforce PRIVATE by default in several modules, expose PUBLIC as needed
move eigen.f90 to obsolete/
(lint) aslong the way
  • 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 5119 2024-07-24 16:46:45Z abarral $
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.