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

Last change on this file since 524 was 512, checked in by acolaitis, 13 years ago

Updated thermals version, that produces quite good overshoots in the inversion region. Series of zmax matches very well the LES and underline one of the weaknesses of convadj. This version will change, as this inversion is only produced in 222 levels. Some resolution tricks must be found for it to work at 32 levels.
File size: 15.9 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)
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      REAL zdz(ngridmx,nlayermx)
68      LOGICAL qtransport_thermals,dtke_thermals
69      INTEGER l,ig,iq,ii(1)
70      CHARACTER (LEN=20) :: modname
71
72!--------------------------------------------------------
73! Local variables for sub-timestep
74!--------------------------------------------------------
75
76      REAL d_t_the(ngridmx,nlayermx), d_q_the(ngridmx,nlayermx,nqmx)
77      REAL d_u_the(ngridmx,nlayermx),d_v_the(ngridmx,nlayermx)
78      REAL dq2_the(ngridmx,nlayermx)
79      INTEGER isplit,nsplit_thermals
80      REAL 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
91!--------------------------------------------------------
92! Diagnostics
93!--------------------------------------------------------
94
95      REAL heatFlux(ngridmx,nlayermx)
96      REAL heatFlux_down(ngridmx,nlayermx)
97      REAL buoyancyOut(ngridmx,nlayermx)
98      REAL buoyancyEst(ngridmx,nlayermx)
99      REAL hfmax(ngridmx),wmax(ngridmx)
100      REAL pbl_teta(ngridmx),dteta(ngridmx,nlayermx)
101
102!--------------------------------------------------------
103! Theta_m
104!--------------------------------------------------------
105
106      INTEGER ico2
107      SAVE ico2
108
109! **********************************************************************
110! Initialization
111! **********************************************************************
112
113      lmax(:)=0.
114      pdu_th(:,:)=0.
115      pdv_th(:,:)=0.
116      pdt_th(:,:)=0.
117      entr_therm(:,:)=0.
118      detr_therm(:,:)=0.
119      q2_therm(:,:)=0.
120      dq2_therm(:,:)=0.
121      ztla(:,:)=0.
122      pbl_dtke(:,:)=0.
123      fm_therm(:,:)=0.
124      zw2(:,:)=0.
125      fraca(:,:)=0.
126      zfraca(:,:)=0.
127      if (tracer) then
128         pdq_th(:,:,:)=0.
129      end if
130      d_t_ajs(:,:)=0.
131      d_u_ajs(:,:)=0.
132      d_v_ajs(:,:)=0.
133      d_q_ajs(:,:,:)=0.
134      heatFlux(:,:)=0.
135      heatFlux_down(:,:)=0.
136      buoyancyOut(:,:)=0.
137      buoyancyEst(:,:)=0.
138      zmaxth(:)=0.
139      lmax_real(:)=0.
140
141
142! **********************************************************************
143! Preparing inputs for the thermals
144! **********************************************************************
145
146       zu(:,:)=pu(:,:)+pdu(:,:)*ptimestep
147       zv(:,:)=pv(:,:)+pdv(:,:)*ptimestep
148       zt(:,:)=pt(:,:)+pdt(:,:)*ptimestep
149
150       pq_therm(:,:,:)=0.
151       qtransport_thermals=.true. !! default setting
152       !call getin("qtransport_thermals",qtransport_thermals)
153
154       if(qtransport_thermals) then
155          if(tracer) then
156                pq_therm(:,:,:)=pq(:,:,:)+pdq(:,:,:)*ptimestep
157          endif
158       endif
159
160!       dtke_thermals=.false. !! default setting
161!       !call getin("dtke_thermals",dtke_thermals)
162!
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     ! ultimately conrols the amount of mass going through the thermals
195                                ! decreasing it increases the thermals effect. Tests at gcm resolution
196                                ! shows that too low values destabilize the model
197                                ! when changing this value, one should check that the surface layer model
198                                ! outputs the correct Cd*u and Ch*u through changing the gustiness coefficient beta
199
200#ifdef MESOSCALE
201         !! valid for timesteps < 200s
202         nsplit_thermals=2
203         r_aspect_thermals=0.7
204#else
205         nsplit_thermals=35
206         r_aspect_thermals=1.5
207#endif
208         call getin("nsplit_thermals",nsplit_thermals)
209         call getin("r_aspect_thermals",r_aspect_thermals)
210
211! **********************************************************************
212! SUB-TIMESTEP LOOP
213! **********************************************************************
214
215         zdt=ptimestep/REAL(nsplit_thermals)
216
217         DO isplit=1,nsplit_thermals
218
219! Initialization of intermediary variables
220
221         zfm_therm(:,:)=0.
222         zentr_therm(:,:)=0.
223         zdetr_therm(:,:)=0.
224         zheatFlux(:,:)=0.
225         zheatFlux_down(:,:)=0.
226         zbuoyancyOut(:,:)=0.
227         zbuoyancyEst(:,:)=0.
228         zzw2(:,:)=0.
229         zmax(:)=0.
230         lmax(:)=0.
231         d_t_the(:,:)=0.
232         d_u_the(:,:)=0.
233         d_v_the(:,:)=0.
234         dq2_the(:,:)=0.
235         if (nqmx .ne. 0) then
236            d_q_the(:,:,:)=0.
237         endif
238
239             CALL thermcell_main_mars(zdt  &
240     &      ,pplay,pplev,pphi,zzlev,zzlay  &
241     &      ,zu,zv,zt,pq_therm,q2_therm  &
242     &      ,d_u_the,d_v_the,d_t_the,d_q_the,dq2_the  &
243     &      ,zfm_therm,zentr_therm,zdetr_therm,lmax,zmax  &
244     &      ,r_aspect_thermals &
245     &      ,zzw2,fraca,zpopsk &
246     &      ,ztla,zheatFlux,zheatFlux_down &
247     &      ,zbuoyancyOut,zbuoyancyEst)
248
249      fact=1./REAL(nsplit_thermals)
250
251            d_t_the(:,:)=d_t_the(:,:)*ptimestep*fact
252!            d_u_the(:,:)=d_u_the(:,:)*fact
253!            d_v_the(:,:)=d_v_the(:,:)*fact
254!            dq2_the(:,:)=dq2_the(:,:)*fact
255            if (ico2 .ne. 0) then
256               d_q_the(:,:,ico2)=d_q_the(:,:,ico2)*fact
257            endif
258
259             zmaxth(:)=zmaxth(:)+zmax(:)*fact
260             lmax_real(:)=lmax_real(:)+float(lmax(:))*fact
261            fm_therm(:,:)=fm_therm(:,:)  &
262     &      +zfm_therm(:,:)*fact
263            entr_therm(:,:)=entr_therm(:,:)  &
264     &       +zentr_therm(:,:)*fact
265            detr_therm(:,:)=detr_therm(:,:)  &
266     &       +zdetr_therm(:,:)*fact
267            zfraca(:,:)=zfraca(:,:) + fraca(:,:)*fact
268
269            heatFlux(:,:)=heatFlux(:,:) &
270     &       +zheatFlux(:,:)*fact
271            heatFlux_down(:,:)=heatFlux_down(:,:) &
272     &       +zheatFlux_down(:,:)*fact
273            buoyancyOut(:,:)=buoyancyOut(:,:) &
274     &       +zbuoyancyOut(:,:)*fact
275            buoyancyEst(:,:)=buoyancyEst(:,:) &
276     &       +zbuoyancyEst(:,:)*fact
277 
278
279            zw2(:,:)=zw2(:,:) + zzw2(:,:)*fact
280
281!  accumulation de la tendance
282
283            d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:)
284!           d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
285!           d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
286            if (ico2 .ne. 0) then
287               d_q_ajs(:,:,ico2)=d_q_ajs(:,:,ico2)+d_q_the(:,:,ico2)
288            endif
289!            dq2_therm(:,:)=dq2_therm(:,:)+dq2_the(:,:)
290!  incrementation des variables meteo
291
292            zt(:,:) = zt(:,:) + d_t_the(:,:)
293!            zu(:,:) = zu(:,:) + d_u_the(:,:)
294!            zv(:,:) = zv(:,:) + d_v_the(:,:)
295            if (ico2 .ne. 0) then
296             pq_therm(:,:,ico2) = &
297     &          pq_therm(:,:,ico2) + d_q_the(:,:,ico2)*ptimestep
298            endif
299!            q2_therm(:,:) = q2_therm(:,:) + dq2_therm(:,:)
300
301
302         ENDDO ! isplit
303!****************************************************************
304
305! Now that we have computed total entrainment and detrainment, we can
306! advect u, v, and q in thermals. (theta already advected). We can do
307! that separatly because u,v,and q are not used in thermcell_main for
308! any thermals-related computation : they are purely passive.
309
310! mass of cells
311      do l=1,nlayermx
312         masse(:,l)=(pplev(:,l)-pplev(:,l+1))/g
313      enddo
314
315! thickness of layers
316      do l=1,nlayermx
317         zdz(:,l)=zzlev(:,l+1)-zzlev(:,l)
318      enddo
319
320      modname='momentum'
321      call thermcell_dqup(ngridmx,nlayermx,ptimestep                &
322     &      ,fm_therm,entr_therm,detr_therm,  &
323     &     masse,zu,d_u_ajs,modname,zdz)
324
325      call thermcell_dqup(ngridmx,nlayermx,ptimestep    &
326     &       ,fm_therm,entr_therm,detr_therm,  &
327     &     masse,zv,d_v_ajs,modname,zdz)
328
329      if (nqmx .ne. 0.) then
330      modname='tracer'
331      DO iq=1,nqmx
332      if (iq .ne. ico2) then
333      call thermcell_dqup(ngridmx,nlayermx,ptimestep     &
334     &     ,fm_therm,entr_therm,detr_therm,  &
335     &    masse,pq_therm(:,:,iq),d_q_ajs(:,:,iq),modname,zdz)
336      endif
337      ENDDO
338      endif
339
340      DO ig=1,ngridmx
341         hfmax(ig)=MAXVAL(heatFlux(ig,:)+heatFlux_down(ig,:))
342         wmax(ig)=MAXVAL(zw2(ig,:))
343      ENDDO
344
345      lmax(:)=nint(lmax_real(:))
346
347! **********************************************************************
348! **********************************************************************
349! **********************************************************************
350! CALLTHERM END
351! **********************************************************************
352! **********************************************************************
353! **********************************************************************
354
355
356! **********************************************************************
357! Preparing outputs
358! **********************************************************************
359
360! Winds and tracers PDU, PDV, and PDQ are DERIVATIVES
361
362           pdu_th(:,:)=d_u_ajs(:,:)
363           pdv_th(:,:)=d_v_ajs(:,:)
364
365           if(qtransport_thermals) then
366              if(tracer) then
367                  pdq_th(:,:,:)=d_q_ajs(:,:,:)
368              endif
369           endif
370
371!           IF(dtke_thermals) THEN
372!              DO l=2,nlayermx
373!                 pbl_dtke(:,l)=0.5*(dq2_therm(:,l-1)+dq2_therm(:,l))
374!              ENDDO
375! 
376!              pbl_dtke(:,1)=0.5*dq2_therm(:,1)
377!              pbl_dtke(:,nlayermx+1)=0.
378!           ENDIF
379
380
381! Temperature PDT is a TENDANCY
382           pdt_th(:,:)=d_t_ajs(:,:)/ptimestep
383
384
385! **********************************************************************
386! Compute the free convection velocity scale for vdifc
387! **********************************************************************
388
389
390! Potential temperature gradient
391
392      dteta(:,nlayermx)=0.
393      DO l=1,nlayermx-1
394         DO ig=1, ngridmx
395            dteta(ig,l) = ((zt(ig,l+1)-zt(ig,l))/zpopsk(ig,l))          &
396     &              /(zzlay(ig,l+1)-zzlay(ig,l))
397         ENDDO
398      ENDDO
399
400! Computation of the pbl mixed layer temperature
401
402      DO ig=1, ngridmx
403         ii=MINLOC(abs(dteta(ig,1:lmax(ig))))
404         pbl_teta(ig) = zt(ig,ii(1))/zpopsk(ig,ii(1))
405      ENDDO
406
407! We follow Spiga et. al 2010 (QJRMS)
408! ------------
409
410      DO ig=1, ngridmx
411         IF (zmax(ig) .gt. 0.) THEN
412            wstar(ig)=(g*zmaxth(ig)*hfmax(ig)/pbl_teta(ig))**(1./3.)
413         ELSE
414            wstar(ig)=0.
415         ENDIF
416      ENDDO
417
418
419
420! **********************************************************************
421! Diagnostics
422! **********************************************************************
423       
424        if(outptherm) then
425        if (ngridmx .eq. 1) then
426        call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',&
427     &                       'kg/m-2',1,entr_therm)
428        call WRITEDIAGFI(ngridmx,'detr_therm','detrainement thermique',&
429     &                       'kg/m-2',1,detr_therm)
430        call WRITEDIAGFI(ngridmx,'fm_therm','flux masse thermique',&
431     &                       'kg/m-2',1,fm_therm)
432        call WRITEDIAGFI(ngridmx,'zw2','vitesse verticale thermique',&
433     &                       'm/s',1,zw2)
434        call WRITEDIAGFI(ngridmx,'heatFlux_up','heatFlux_updraft',&
435     &                       'SI',1,heatFlux)
436       call WRITEDIAGFI(ngridmx,'heatFlux_down','heatFlux_downdraft',&
437     &                       'SI',1,heatFlux_down)
438        call WRITEDIAGFI(ngridmx,'fraca','fraction coverage',&
439     &                       'percent',1,fraca)
440        call WRITEDIAGFI(ngridmx,'buoyancyOut','buoyancyOut',&
441     &                       'm.s-2',1,buoyancyOut)
442        call WRITEDIAGFI(ngridmx,'buoyancyEst','buoyancyEst',&
443     &                       'm.s-2',1,buoyancyEst)
444        call WRITEDIAGFI(ngridmx,'d_t_th',  &
445     &         'tendance temp TH','K',1,d_t_ajs)
446        call WRITEDIAGFI(ngridmx,'zmax',  &
447     &         'pbl height','m',0,zmaxth)
448      else
449
450        call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',&
451     &                       'kg/m-2',3,entr_therm)
452        call WRITEDIAGFI(ngridmx,'detr_therm','detrainement thermique',&
453     &                       'kg/m-2',3,detr_therm)
454        call WRITEDIAGFI(ngridmx,'fm_therm','flux masse thermique',&
455     &                       'kg/m-2',3,fm_therm)
456        call WRITEDIAGFI(ngridmx,'zw2','vitesse verticale thermique',&
457     &                       'm/s',3,zw2)
458        call WRITEDIAGFI(ngridmx,'heatFlux','heatFlux',&
459     &                       'SI',3,heatFlux)
460        call WRITEDIAGFI(ngridmx,'buoyancyOut','buoyancyOut',&
461     &                       'SI',3,buoyancyOut)
462        call WRITEDIAGFI(ngridmx,'d_t_th',  &
463     &         'tendance temp TH','K',3,d_t_ajs)
464
465      endif
466      endif
467
468       END
Note: See TracBrowser for help on using the repository browser.