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

Last change on this file since 508 was 508, checked in by acolaitis, 14 years ago

Added molar-mass vertical gradient computations in thermals, coupled to the advection scheme for co2 tracer of thermals. As a result, there are now gentle thermals in the polar night (associated with a wstar of about 1.5 m.s-1).

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