source: LMDZ6/branches/Ocean_skin/libf/phylmd/thermcell_flux2.F90 @ 5185

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