source: trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90 @ 640

Last change on this file since 640 was 628, checked in by acolaitis, 13 years ago

Correction to advection problems in thermals + made the thermals model faster by limiting the vertical extension of loops with the height reached by thermals.

File size: 17.8 KB
Line 
1!
2! AC 2011-01-05
3!
4      SUBROUTINE calltherm_interface (firstcall, &
5     & long,lati,zzlev,zzlay, &
6     & ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2, &
7     & pplay,pplev,pphi,zpopsk, &
8     & pdu_th,pdv_th,pdt_th,pdq_th,lmax,zmaxth,pbl_dtke,hfmax,wstar)
9
10       USE ioipsl_getincom
11
12      implicit none
13#include "callkeys.h"
14#include "dimensions.h"
15#include "dimphys.h"
16#include "comcstfi.h"
17#include "tracer.h"
18
19!--------------------------------------------------------
20! Input Variables
21!--------------------------------------------------------
22
23      REAL, INTENT(IN) :: ptimestep
24      REAL, INTENT(IN) :: pplev(ngridmx,nlayermx+1),pplay(ngridmx,nlayermx)
25      REAL, INTENT(IN) :: pphi(ngridmx,nlayermx)
26      REAL, INTENT(IN) :: pu(ngridmx,nlayermx),pv(ngridmx,nlayermx)
27      REAL, INTENT(IN) :: pt(ngridmx,nlayermx),pq(ngridmx,nlayermx,nqmx)
28      REAL, INTENT(IN) :: zzlay(ngridmx,nlayermx)
29      REAL, INTENT(IN) :: zzlev(ngridmx,nlayermx+1)
30      LOGICAL, INTENT(IN) :: firstcall
31      REAL, INTENT(IN) :: pdu(ngridmx,nlayermx),pdv(ngridmx,nlayermx)
32      REAL, INTENT(IN) :: pdq(ngridmx,nlayermx,nqmx),pdt(ngridmx,nlayermx)
33      REAL, INTENT(IN) :: q2(ngridmx,nlayermx+1)
34      REAL, INTENT(IN) :: long(ngridmx),lati(ngridmx)
35      REAL, INTENT(IN) :: zpopsk(ngridmx,nlayermx)
36
37!--------------------------------------------------------
38! Output Variables
39!--------------------------------------------------------
40
41      REAL, INTENT(OUT) :: pdu_th(ngridmx,nlayermx)
42      REAL, INTENT(OUT) :: pdv_th(ngridmx,nlayermx)
43      REAL, INTENT(OUT) :: pdt_th(ngridmx,nlayermx)
44      REAL, INTENT(OUT) :: pdq_th(ngridmx,nlayermx,nqmx)
45      INTEGER, INTENT(OUT) :: lmax(ngridmx)
46      REAL, INTENT(OUT) :: zmaxth(ngridmx)
47      REAL, INTENT(OUT) :: pbl_dtke(ngridmx,nlayermx+1)
48      REAL, INTENT(OUT) :: wstar(ngridmx)
49
50!--------------------------------------------------------
51! Thermals local variables
52!--------------------------------------------------------
53      REAL zu(ngridmx,nlayermx), zv(ngridmx,nlayermx)
54      REAL zt(ngridmx,nlayermx)
55      REAL d_t_ajs(ngridmx,nlayermx)
56      REAL d_u_ajs(ngridmx,nlayermx), d_q_ajs(ngridmx,nlayermx,nqmx)
57      REAL d_v_ajs(ngridmx,nlayermx)
58      REAL fm_therm(ngridmx,nlayermx+1), entr_therm(ngridmx,nlayermx)
59      REAL detr_therm(ngridmx,nlayermx),detrmod(ngridmx,nlayermx)
60      REAL zw2(ngridmx,nlayermx+1)
61      REAL fraca(ngridmx,nlayermx+1),zfraca(ngridmx,nlayermx+1)
62      REAL ztla(ngridmx,nlayermx)
63      REAL q_therm(ngridmx,nlayermx), pq_therm(ngridmx,nlayermx,nqmx)
64      REAL q2_therm(ngridmx,nlayermx), dq2_therm(ngridmx,nlayermx)
65      REAL lmax_real(ngridmx)
66      REAL masse(ngridmx,nlayermx)
67      LOGICAL qtransport_thermals,dtke_thermals
68      INTEGER l,ig,iq,ii(1)
69      CHARACTER (LEN=20) modname
70
71!--------------------------------------------------------
72! Local variables for sub-timestep
73!--------------------------------------------------------
74
75      REAL d_t_the(ngridmx,nlayermx), d_q_the(ngridmx,nlayermx,nqmx)
76      REAL d_u_the(ngridmx,nlayermx),d_v_the(ngridmx,nlayermx)
77      REAL dq2_the(ngridmx,nlayermx)
78      INTEGER isplit
79      INTEGER,SAVE :: nsplit_thermals
80      REAL, SAVE :: r_aspect_thermals
81      REAL fact
82      REAL zfm_therm(ngridmx,nlayermx+1),zdt
83      REAL zentr_therm(ngridmx,nlayermx),zdetr_therm(ngridmx,nlayermx)
84      REAL zheatFlux(ngridmx,nlayermx)
85      REAL zheatFlux_down(ngridmx,nlayermx)
86      REAL zbuoyancyOut(ngridmx,nlayermx)
87      REAL zbuoyancyEst(ngridmx,nlayermx)
88      REAL zzw2(ngridmx,nlayermx+1)
89      REAL zmax(ngridmx)
90      INTEGER ndt,zlmax
91
92!--------------------------------------------------------
93! Diagnostics
94!--------------------------------------------------------
95
96      REAL heatFlux(ngridmx,nlayermx)
97      REAL heatFlux_down(ngridmx,nlayermx)
98      REAL buoyancyOut(ngridmx,nlayermx)
99      REAL buoyancyEst(ngridmx,nlayermx)
100      REAL hfmax(ngridmx),wmax(ngridmx)
101      REAL pbl_teta(ngridmx),dteta(ngridmx,nlayermx)
102
103!--------------------------------------------------------
104! Theta_m
105!--------------------------------------------------------
106
107      INTEGER ico2
108      SAVE ico2
109
110! **********************************************************************
111! Initialization
112! **********************************************************************
113
114      lmax(:)=0
115      pdu_th(:,:)=0.
116      pdv_th(:,:)=0.
117      pdt_th(:,:)=0.
118      entr_therm(:,:)=0.
119      detr_therm(:,:)=0.
120      q2_therm(:,:)=0.
121      dq2_therm(:,:)=0.
122      ztla(:,:)=0.
123      pbl_dtke(:,:)=0.
124      fm_therm(:,:)=0.
125      zw2(:,:)=0.
126      fraca(:,:)=0.
127      zfraca(:,:)=0.
128      if (tracer) then
129         pdq_th(:,:,:)=0.
130      end if
131      d_t_ajs(:,:)=0.
132      d_u_ajs(:,:)=0.
133      d_v_ajs(:,:)=0.
134      d_q_ajs(:,:,:)=0.
135      heatFlux(:,:)=0.
136      heatFlux_down(:,:)=0.
137      buoyancyOut(:,:)=0.
138      buoyancyEst(:,:)=0.
139      zmaxth(:)=0.
140      lmax_real(:)=0.
141
142
143! **********************************************************************
144! Preparing inputs for the thermals
145! **********************************************************************
146
147       zu(:,:)=pu(:,:)+pdu(:,:)*ptimestep
148       zv(:,:)=pv(:,:)+pdv(:,:)*ptimestep
149       zt(:,:)=pt(:,:)+pdt(:,:)*ptimestep
150
151       pq_therm(:,:,:)=0.
152       qtransport_thermals=.true. !! default setting
153       !call getin("qtransport_thermals",qtransport_thermals)
154
155       if(qtransport_thermals) then
156          if(tracer) then
157                pq_therm(:,:,:)=pq(:,:,:)+pdq(:,:,:)*ptimestep
158          endif
159       endif
160
161       dtke_thermals=.false. !! default setting
162       call getin("dtke_thermals",dtke_thermals)
163       IF(dtke_thermals) THEN
164          DO l=1,nlayermx
165              q2_therm(:,l)=0.5*(q2(:,l)+q2(:,l+1))
166          ENDDO
167       ENDIF
168
169! **********************************************************************
170! Polar night mixing : theta_m
171! **********************************************************************
172
173      if(firstcall) then
174        ico2=0
175        if (tracer) then
176!     Prepare Special treatment if one of the tracers is CO2 gas
177           do iq=1,nqmx
178             if (noms(iq).eq."co2") then
179                ico2=iq
180             end if
181           enddo
182        endif
183      endif !of if firstcall
184
185
186! **********************************************************************
187! **********************************************************************
188! **********************************************************************
189! CALLTHERM
190! **********************************************************************
191! **********************************************************************
192! **********************************************************************
193
194!         r_aspect_thermals     ! Mainly control the shape of the temperature profile
195                                ! in the surface layer. Decreasing it goes toward
196                                ! a convective-adjustment like profile.
197!         nsplit_thermals       ! Sub-timestep for the thermals. Very dependant on the
198                                ! chosen timestep for the radiative transfer.
199                                ! It is recommended to run with 96 timestep per day and
200                                ! iradia = 1., configuration in which thermals can run
201                                ! very well with a sub-timestep of 10.
202         IF (firstcall) THEN
203#ifdef MESOSCALE
204            !! valid for timesteps < 200s
205            nsplit_thermals=4
206            r_aspect_thermals=0.7
207#else
208            IF ((ptimestep .le. 3699.*24./96.) .and. (iradia .eq. 1)) THEN
209               nsplit_thermals=10
210            ELSE
211               nsplit_thermals=35
212            ENDIF
213            r_aspect_thermals=1.
214#endif
215            call getin("nsplit_thermals",nsplit_thermals)
216            call getin("r_aspect_thermals",r_aspect_thermals)
217         ENDIF
218
219! **********************************************************************
220! SUB-TIMESTEP LOOP
221! **********************************************************************
222
223         zdt=ptimestep/REAL(nsplit_thermals)
224
225         DO isplit=1,nsplit_thermals
226
227! Initialization of intermediary variables
228
229!         zfm_therm(:,:)=0. !init is done inside
230!         zentr_therm(:,:)=0.
231!         zdetr_therm(:,:)=0.
232!         zheatFlux(:,:)=0. 
233!         zheatFlux_down(:,:)=0.
234!         zbuoyancyOut(:,:)=0.
235!         zbuoyancyEst(:,:)=0.
236         zzw2(:,:)=0.
237         zmax(:)=0.
238         lmax(:)=0
239!         d_t_the(:,:)=0. !init is done inside
240
241!         d_u_the(:,:)=0. !transported outside
242!         d_v_the(:,:)=0.
243         dq2_the(:,:)=0.
244
245         if (nqmx .ne. 0 .and. ico2 .ne. 0) then
246            d_q_the(:,:,ico2)=0.
247         endif
248
249             CALL thermcell_main_mars(zdt  &
250     &      ,pplay,pplev,pphi,zzlev,zzlay  &
251     &      ,zu,zv,zt,pq_therm,q2_therm  &
252     &      ,d_u_the,d_v_the,d_t_the,d_q_the,dq2_the  &
253     &      ,zfm_therm,zentr_therm,zdetr_therm,lmax,zmax  &
254     &      ,r_aspect_thermals &
255     &      ,zzw2,fraca,zpopsk &
256     &      ,ztla,zheatFlux,zheatFlux_down &
257     &      ,zbuoyancyOut,zbuoyancyEst)
258
259      fact=1./REAL(nsplit_thermals)
260
261            d_t_the(:,:)=d_t_the(:,:)*ptimestep*fact
262!            d_u_the(:,:)=d_u_the(:,:)*ptimestep*fact
263!            d_v_the(:,:)=d_v_the(:,:)*ptimestep*fact
264            dq2_the(:,:)=dq2_the(:,:)*fact
265            if (ico2 .ne. 0) then
266               d_q_the(:,:,ico2)=d_q_the(:,:,ico2)*ptimestep*fact
267            endif
268
269            zmaxth(:)=zmaxth(:)+zmax(:)*fact
270            lmax_real(:)=lmax_real(:)+float(lmax(:))*fact
271            fm_therm(:,:)=fm_therm(:,:)  &
272     &      +zfm_therm(:,:)*fact
273            entr_therm(:,:)=entr_therm(:,:)  &
274     &       +zentr_therm(:,:)*fact
275            detr_therm(:,:)=detr_therm(:,:)  &
276     &       +zdetr_therm(:,:)*fact
277            zfraca(:,:)=zfraca(:,:) + fraca(:,:)*fact
278
279            heatFlux(:,:)=heatFlux(:,:) &
280     &       +zheatFlux(:,:)*fact
281            heatFlux_down(:,:)=heatFlux_down(:,:) &
282     &       +zheatFlux_down(:,:)*fact
283            buoyancyOut(:,:)=buoyancyOut(:,:) &
284     &       +zbuoyancyOut(:,:)*fact
285            buoyancyEst(:,:)=buoyancyEst(:,:) &
286     &       +zbuoyancyEst(:,:)*fact
287 
288
289            zw2(:,:)=zw2(:,:) + zzw2(:,:)*fact
290
291!  accumulation de la tendance
292
293           d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:)
294!           d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
295!           d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
296            if (ico2 .ne. 0) then
297               d_q_ajs(:,:,ico2)=d_q_ajs(:,:,ico2)+d_q_the(:,:,ico2)
298            endif
299!            dq2_therm(:,:)=dq2_therm(:,:)+dq2_the(:,:)
300!  incrementation des variables meteo
301
302            zt(:,:) = zt(:,:) + d_t_the(:,:)
303!            zu(:,:) = zu(:,:) + d_u_the(:,:)
304!            zv(:,:) = zv(:,:) + d_v_the(:,:)
305            if (ico2 .ne. 0) then
306             pq_therm(:,:,ico2) = &
307     &          pq_therm(:,:,ico2) + d_q_the(:,:,ico2)
308            endif
309!            q2_therm(:,:) = q2_therm(:,:) + dq2_therm(:,:)
310
311
312         ENDDO ! isplit
313!****************************************************************
314
315      lmax(:)=nint(lmax_real(:))
316      zlmax=MAXVAL(lmax(:))+2
317      if (zlmax .ge. nlayermx) then
318        print*,'thermals have reached last layer of the model'
319        print*,'this is not good !'
320      endif
321
322
323! Now that we have computed total entrainment and detrainment, we can
324! advect u, v, and q in thermals. (theta already advected). We can do
325! that separatly because u,v,and q are not used in thermcell_main for
326! any thermals-related computation : they are purely passive.
327
328! mass of cells
329      do l=1,nlayermx
330         masse(:,l)=(pplev(:,l)-pplev(:,l+1))/g
331      enddo
332
333      detrmod(:,:)=0.
334      do l=1,zlmax
335         do ig=1,ngridmx
336            detrmod(ig,l)=fm_therm(ig,l)-fm_therm(ig,l+1) &
337     &      +entr_therm(ig,l)
338            if (detrmod(ig,l).lt.0.) then
339               entr_therm(ig,l)=entr_therm(ig,l)-detrmod(ig,l)
340               detrmod(ig,l)=0.
341            endif
342         enddo
343      enddo
344      ndt=10
345      call thermcell_dqup(ngridmx,nlayermx,ptimestep                &
346     &      ,fm_therm,entr_therm,detrmod,  &
347     &     masse,zu,d_u_ajs,ndt,zlmax)
348
349      call thermcell_dqup(ngridmx,nlayermx,ptimestep    &
350     &       ,fm_therm,entr_therm,detrmod,  &
351     &     masse,zv,d_v_ajs,ndt,zlmax)
352
353      if (nqmx .ne. 0.) then
354      DO iq=1,nqmx
355      if (iq .ne. ico2) then
356      call thermcell_dqup(ngridmx,nlayermx,ptimestep     &
357     &     ,fm_therm,entr_therm,detrmod,  &
358     &    masse,pq_therm(:,:,iq),d_q_ajs(:,:,iq),ndt,zlmax)
359      endif
360      ENDDO
361      endif
362
363      if (dtke_thermals) then
364      detrmod(:,:)=0.
365      ndt=1
366      do l=1,zlmax
367         do ig=1,ngridmx
368            detrmod(ig,l)=fm_therm(ig,l)-fm_therm(ig,l+1) &
369     &      +entr_therm(ig,l)
370            if (detrmod(ig,l).lt.0.) then
371               entr_therm(ig,l)=entr_therm(ig,l)-detrmod(ig,l)
372               detrmod(ig,l)=0.
373            endif
374         enddo
375      enddo
376      call thermcell_dqup(ngridmx,nlayermx,ptimestep     &
377     &     ,fm_therm,entr_therm,detrmod,  &
378     &    masse,q2_therm,dq2_therm,ndt,zlmax)
379      endif
380
381      DO ig=1,ngridmx
382         hfmax(ig)=MAXVAL(heatFlux(ig,:)+heatFlux_down(ig,:))
383         wmax(ig)=MAXVAL(zw2(ig,:))
384      ENDDO
385
386! **********************************************************************
387! **********************************************************************
388! **********************************************************************
389! CALLTHERM END
390! **********************************************************************
391! **********************************************************************
392! **********************************************************************
393
394
395! **********************************************************************
396! Preparing outputs
397! **********************************************************************
398
399      do l=1,zlmax
400        pdu_th(:,l)=d_u_ajs(:,l)
401        pdv_th(:,l)=d_v_ajs(:,l)
402      enddo
403
404           if(qtransport_thermals) then
405              if(tracer) then
406               do iq=1,nqmx
407                if (iq .ne. ico2) then
408                  do l=1,zlmax
409                     pdq_th(:,l,iq)=d_q_ajs(:,l,iq)
410                  enddo
411                else
412                  do l=1,zlmax
413                     pdq_th(:,l,iq)=d_q_ajs(:,l,iq)/ptimestep
414                  enddo
415                endif
416               enddo
417              endif
418           endif
419
420           IF(dtke_thermals) THEN
421              DO l=2,nlayermx
422                 pbl_dtke(:,l)=0.5*(dq2_therm(:,l-1)+dq2_therm(:,l))
423              ENDDO
424 
425              pbl_dtke(:,1)=0.5*dq2_therm(:,1)
426              pbl_dtke(:,nlayermx+1)=0.
427           ENDIF
428
429           do l=1,zlmax
430              pdt_th(:,l)=d_t_ajs(:,l)/ptimestep
431           enddo
432
433
434! **********************************************************************
435! Compute the free convection velocity scale for vdifc
436! **********************************************************************
437
438
439! Potential temperature gradient
440
441      dteta(:,nlayermx)=0.
442      DO l=1,nlayermx-1
443         DO ig=1, ngridmx
444            dteta(ig,l) = ((zt(ig,l+1)-zt(ig,l))/zpopsk(ig,l))          &
445     &              /(zzlay(ig,l+1)-zzlay(ig,l))
446         ENDDO
447      ENDDO
448
449! Computation of the pbl mixed layer temperature
450
451      DO ig=1, ngridmx
452         ii=MINLOC(abs(dteta(ig,1:lmax(ig))))
453         pbl_teta(ig) = zt(ig,ii(1))/zpopsk(ig,ii(1))
454      ENDDO
455
456! We follow Spiga et. al 2010 (QJRMS)
457! ------------
458
459      DO ig=1, ngridmx
460         IF (zmax(ig) .gt. 0.) THEN
461            wstar(ig)=(g*zmaxth(ig)*hfmax(ig)/pbl_teta(ig))**(1./3.)
462         ELSE
463            wstar(ig)=0.
464         ENDIF
465      ENDDO
466
467
468
469! **********************************************************************
470! Diagnostics
471! **********************************************************************
472       
473        if(outptherm) then
474        if (ngridmx .eq. 1) then
475        call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',&
476     &                       'kg/m-2',1,entr_therm)
477        call WRITEDIAGFI(ngridmx,'detr_therm','detrainement thermique',&
478     &                       'kg/m-2',1,detr_therm)
479        call WRITEDIAGFI(ngridmx,'fm_therm','flux masse thermique',&
480     &                       'kg/m-2',1,fm_therm)
481        call WRITEDIAGFI(ngridmx,'zw2','vitesse verticale thermique',&
482     &                       'm/s',1,zw2)
483        call WRITEDIAGFI(ngridmx,'heatFlux_up','heatFlux_updraft',&
484     &                       'SI',1,heatFlux)
485       call WRITEDIAGFI(ngridmx,'heatFlux_down','heatFlux_downdraft',&
486     &                       'SI',1,heatFlux_down)
487        call WRITEDIAGFI(ngridmx,'fraca','fraction coverage',&
488     &                       'percent',1,fraca)
489        call WRITEDIAGFI(ngridmx,'buoyancyOut','buoyancyOut',&
490     &                       'm.s-2',1,buoyancyOut)
491        call WRITEDIAGFI(ngridmx,'buoyancyEst','buoyancyEst',&
492     &                       'm.s-2',1,buoyancyEst)
493        call WRITEDIAGFI(ngridmx,'d_t_th',  &
494     &         'tendance temp TH','K',1,d_t_ajs)
495        call WRITEDIAGFI(ngridmx,'d_q_th',  &
496     &         'tendance traceur TH','kg/kg',1,d_q_ajs)
497        call WRITEDIAGFI(ngridmx,'zmax',  &
498     &         'pbl height','m',0,zmaxth)
499        call WRITEDIAGFI(ngridmx,'d_u_th',  &
500     &         'tendance moment','m/s',1,pdu_th)
501      else
502
503        call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',&
504     &                       'kg/m-2',3,entr_therm)
505        call WRITEDIAGFI(ngridmx,'detr_therm','detrainement thermique',&
506     &                       'kg/m-2',3,detr_therm)
507        call WRITEDIAGFI(ngridmx,'fm_therm','flux masse thermique',&
508     &                       'kg/m-2',3,fm_therm)
509        call WRITEDIAGFI(ngridmx,'zw2','vitesse verticale thermique',&
510     &                       'm/s',3,zw2)
511        call WRITEDIAGFI(ngridmx,'heatFlux','heatFlux',&
512     &                       'SI',3,heatFlux)
513        call WRITEDIAGFI(ngridmx,'buoyancyOut','buoyancyOut',&
514     &                       'SI',3,buoyancyOut)
515        call WRITEDIAGFI(ngridmx,'d_t_th',  &
516     &         'tendance temp TH','K',3,d_t_ajs)
517
518      endif
519      endif
520
521       END
Note: See TracBrowser for help on using the repository browser.