source: LMDZ6/branches/DYNAMICO-conv-GC/libf/phymar/PHY_Atm_CP_RUN.f90 @ 4249

Last change on this file since 4249 was 2089, checked in by Laurent Fairhead, 10 years ago

Inclusion de la physique de MAR


Integration of MAR physics

File size: 18.6 KB
Line 
1      subroutine PHY_Atm_CP_RUN(mzc,kcolc)
2
3!------------------------------------------------------------------------------+
4!                                                         Mon 17-Jun-2013  MAR |
5!     subroutine PHY_Atm_CP_RUN    interfaces                                  |
6!                Bechtold et al. (2001) Convective Parameterization with MAR   |
7!                                                                              |
8!     version 3.p.4.1 created by H. Gallee,               Mon  8-Apr-2013      |
9!           Last Modification by H. Gallee,               Mon 17-Jun-2013      |
10!                                                                              |
11!------------------------------------------------------------------------------+
12!                                                                              |
13!     INPUT                                                                    |
14!     ^^^^^        it_EXP             : Experiment Iteration Counter           |
15!                  dt__CP             : Mass Flux  Scheme:   Time Step         |
16!                  dxHOST             : grid spacing  of  HOST MODEL       [m] |
17!                  Ta__DY(kcolp,mzpp) : air  temperature                   [K] |
18!                  qs__CM(kcolp,mzpp) : air  snow  Particl. concentr.  [kg/kg] |
19!                  qr__CM(kcolp,mzp ) : air  rain  drops    concentr.  [kg/kg] |
20!                                                                              |
21!     INPUT/OUTPUT                                                             |
22!     ^^^^^^^^^^^^ qv__DY(kcolp,mzpp) : Specific Humidity              [kg/kg] |
23!                  qw__CM(kcolp,mzp ) : air  cloud droplets concentr.  [kg/kg] |
24!                  qi__CM(kcolp,mzp ) : air  cloud crystals concentr.  [kg/kg] |
25!                  snowCP(kcolp     ) : Snow (convective)                  [m] |
26!                  snowCM(kcolp     ) : Snow (convective + stratiform)     [m] |
27!                  rainCP(kcolp     ) : Rain (convective)                  [m] |
28!                  rainCM(kcolp     ) : Rain (convective + stratiform)     [m] |
29!                                                                              |
30!     OUTPUT       dpktCP(kcolp,mzp ) : Reduc. Pot.Temperat.Tendency   [K/X/s] |
31!     ^^^^^^       dqv_CP(kcolp,mzp ) : Specific   Humidity Tendency [kg/kg/s] |
32!                  dqw_CP(kcolp,mzp ) : cloud dropl.Concent.Tendency [kg/kg/s] |
33!                  dqi_CP(kcolp,mzp ) : cloud cryst.Concent.Tendency [kg/kg/s] |
34!                                                                              |
35!                  dss_CP(kcolp     ) : Snow (convective)   Tendency     [m/s] |
36!                  drr_CP(kcolp     ) : Rain (convective)   Tendency     [m/s] |
37!                                                                              |
38!                  pkt_DY(kcolp,mzp ) : Reduced  Potential Temperature   [K/X] |
39!                                                                              |
40!                  CAPECP(kcolp     ) : Convective Avail.Potent.Energy         |
41!                                                                              |
42!     REFER. :  MesoNH MASS FLUX CONVECTION Routine                            |
43!     ^^^^^^^^         (Bechtold et al., 2001, QJRMS 127, pp 869-886)          |
44!                                                                              |
45!   # OPTION : #EW  Energy and Water  Conservation                             |
46!   # ^^^^^^^^                                                                 |
47!                                                                              |
48!     MODIF. HGallee: 18-11-2004: Adaptation to CVAmnh.f90.laurent             |
49!     ^^^^^^                      (Argument kensbl of CONVECTION removed)      |
50!                                                                              |
51!------------------------------------------------------------------------------+
52
53      use Mod_Real
54      use Mod_PHY____dat
55      use Mod_PHY____grd
56      use Mod_PHY_CP_ctr
57      use Mod_PHY_CP_dat
58      use Mod_PHY_CP_grd
59      use Mod_PHY_CP_kkl
60      use Mod_PHY_DY_kkl
61      use Mod_PHY_AT_kkl
62      use Mod_PHY_CM_kkl
63
64      use Mod_Atm_CP_RUN
65
66      IMPLICIT NONE
67
68
69
70! Global Variables
71! ================
72
73      integer                                       ::  mzc                      !  Nb of levels
74      integer                                       ::  kcolc                    !  Nb of columns
75
76
77
78
79! Local  Variables
80! ================
81
82      real(kind=real8)                              ::  bANA                     !
83      real(kind=real8)                              ::  zANA                     !
84
85      integer                                       ::  i                        !  x-Axis Index
86      integer                                       ::  j                        !  y-Axis Index
87      integer                                       ::  k                        !  Level  Index (from top    to bottom)
88      integer                                       ::  klc                      !  Level  Index (from bottom to top   )
89      integer                                       ::  ikl                      !  Column Index
90
91      REAL                                          ::  Pdxdy0(kcolc)
92      REAL                                          ::  P_pa_0(kcolc,mzc)
93      REAL                                          ::  P_za_0(kcolc,mzc)
94      REAL                                          ::  P_Ta_0(kcolc,mzc)
95      REAL                                          ::  P_Qa_0(kcolc,mzc)
96      REAL                                          ::  P_Qw_0(kcolc,mzc)
97      REAL                                          ::  P_Qi_0(kcolc,mzc)
98      REAL                                          ::  P_Ua_0(kcolc,mzc)
99      REAL                                          ::  P_Va_0(kcolc,mzc)
100      REAL                                          ::  P_Wa_0(kcolc,mzc)
101
102      integer                                       ::  Kstep1(kcolc)            ! convective counter
103      integer                                       ::  K_CbT1(kcolc)            ! cloud top  level
104      integer                                       ::  K_CbB1(kcolc)            ! cloud base level
105
106      integer, parameter                            ::  KTCCH0=1                 !
107      REAL                                          ::  P_CH_0(kcolc,mzc,KTCCH0)
108      REAL                                          ::  PdCH_1(kcolc,mzc,KTCCH0)
109
110      REAL                                          ::  PdTa_1(kcolc,mzc)
111      REAL                                          ::  PdQa_1(kcolc,mzc)
112      REAL                                          ::  PdQw_1(kcolc,mzc)
113      REAL                                          ::  PdQi_1(kcolc,mzc)
114      REAL                                          ::  Pdrr_1(kcolc)
115      REAL                                          ::  Pdss_1(kcolc)
116      REAL                                          ::  PuMF_1(kcolc,mzc)        !   Upward        Mass Flux
117      REAL                                          ::  PdMF_1(kcolc,mzc)        ! Downward        Mass Flux
118      REAL                                          ::  Pfrr_1(kcolc,mzc)        ! Liquid Precipitation Flux
119      REAL                                          ::  Pfss_1(kcolc,mzc)        !  Solid Precipitation Flux
120      REAL                                          ::  Pcape1(kcolc)            !   CAPE [J/kg]
121
122
123!  Diagnostic Variables
124!  --------------------
125
126! #EW integer                                       ::  irmx  ,jrmx  ,iter_0     !
127! #EW real(kind=real8)                              ::  rr_max,temp_r,energ0     !
128! #EW real(kind=real8)                              ::  water0,waterb            !
129
130
131
132!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
133!                                                                                !
134!  ALLOCATION
135!  ----------
136
137      IF (it_RUN.EQ.1 .OR. FlagDALLOC)                                      THEN !
138          allocate  ( wa_ANA(kcolc,mzc) )                                        !  ANAbatic       Wind        Speed     [m/s]
139      END IF
140!                                                                                !
141!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
142
143
144
145
146! Update Convective Mass Flux
147! ===========================
148
149      IF     ( mod(iitCV0,jjtCV0).EQ.0 )                            THEN
150
151! Martin control
152        PRINT*,'Dans PHY_Atm_CP_RUN'
153        PRINT*,'size(ii__AP)',size(ii__AP)
154        PRINT*,'size(jj__AP)',size(jj__AP)
155        PRINT*,'size(wa_ANA)=',size(wa_ANA)
156        PRINT*,'ii__AP(1)=',ii__AP(1)
157        PRINT*,'jj__AP(1)=',jj__AP(1)
158        PRINT*,'ii__AP(kcolp)=',ii__AP(kcolp)
159        PRINT*,'jj__AP(kcolp)=',jj__AP(kcolp)
160        PRINT*,'mzp=',mzp
161! Martin control
162
163        DO       ikl = 1,kcolp
164
165!           PRINT*,'ikl=',ikl
166
167           i   = ii__AP(ikl)
168           j   = jj__AP(ikl)
169
170!          PRINT*,'ii__AP(',ikl,')=',ii__AP(ikl)
171!          PRINT*,'jj__AP(',ikl,')=',jj__AP(ikl)
172
173
174!  Contribution from Subgrid Mountain Breeze
175!  -----------------------------------------
176
177            DO k=1,mzp
178                  wa_ANA(ikl,k   ) =  0.0
179            END DO
180          IF(Lo_ANA)                                                THEN
181            DO k=1,mzp
182                  bANA             =  min(Z___DY(ikl,k  ),  zi__AT(ikl))
183                  zANA             =        hANA(ikl)      +  bANA * 2.0
184              IF (Z___DY(ikl,k   )   .LE.   zANA             .AND.      &
185     &            Ta__DY(ikl,mzpp)   .GT. Ta__DY(ikl,mzp))          THEN
186                  wa_ANA(ikl,k   ) =        rANA           *  bANA * 0.5 !  Half Integrated Horizontal Divergence
187
188              END IF
189            END DO
190          END IF
191
192
193
194!  Mass Flux convective Scheme: Set Up Vertical Profiles
195!  -----------------------------------------------------
196
197             Pdxdy0(ikl)      = dxHOST * dxHOST                          !  grid area                           [m2]
198          DO klc= 1,mzp
199             k =    mzpp-klc
200             P_pa_0(ikl,klc) = (psa_DY(  ikl)*sigma(k)  + pt__DY) *1.e3  !  pressure     in layer               [Pa]
201             P_za_0(ikl,klc) =  Z___DY(ikl,k)                            !  height of model layer                [m]
202             P_Ta_0(ikl,klc) =  Ta__DY(ikl,k)                            !  grid scale T           at time t     [K]
203             P_Qa_0(ikl,klc) =  qv__DY(ikl,k)                            !  grid scale water vapor at time t [kg/kg]
204             P_Qw_0(ikl,klc) =  qw__CM(ikl,k) / (1.0-qw__CM(ikl,k))      !  grid scale Cloud drops at time t [kg/kg]
205             P_Qi_0(ikl,klc) =  qi__CM(ikl,k) / (1.0-qi__CM(ikl,k))      !  grid scale Cloud ice   at time t [kg/kg]
206             P_Ua_0(ikl,klc) =  ua__DY(ikl,k)                            !  grid scale hor. wind u at time t   [m/s]
207             P_Va_0(ikl,klc) =  va__DY(ikl,k)                            !  grid scale hor. wind v at time t   [m/s]
208             P_Wa_0(ikl,klc) =  wa__DY(ikl,k) +      wa_ANA(ikl,k)       !  grid scale vertic.wind at time t   [m/s]
209          END DO
210
211        END DO ! ikl = 1,kcolp
212
213
214!  Mass Flux convective Scheme: Bechtold et al. (2001) Convective Parameterization
215!  -------------------------------------------------------------------------------
216
217!         ***************
218          call CONVECTION(                                               &
219     &             kcolc , mzc   , kidia0, kfdia0, kbdia0, ktdia0,       &
220     &             pdtCV , Odeep , Oshal , Orset0, Odown0, kIce_0,       &
221     &             OsetA0, PTdcv , PTscv ,                               &
222     &             kensbl,                                               &
223     &     P_pa_0, P_za_0, Pdxdy0,                                       &
224     &     P_Ta_0, P_Qa_0, P_Qw_0, P_Qi_0, P_Ua_0, P_Va_0, P_Wa_0,       &
225     &     Kstep1, PdTa_1, PdQa_1, PdQw_1, PdQi_1,                       &
226     &                     Pdrr_1, Pdss_1,                               &
227     &     PuMF_1, PdMF_1, Pfrr_1, Pfss_1, Pcape1, K_CbT1, K_CbB1,       &
228     &     OCvTC0, KTCCH0, P_CH_0, PdCH_1)
229!         ***************
230
231
232
233!  Mass Flux convective Scheme: products
234!  -------------------------------------
235
236        DO       ikl = 1,kcolp
237             CAPECP(ikl)   = Pcape1(ikl)
238             timeCP(ikl)   = Kstep1(ikl)     *dt__CP
239             drr_CP(ikl)   = Pdrr_1(ikl)
240             dss_CP(ikl)   = Pdss_1(ikl)
241
242          DO klc= 1,mzp
243             k  =   mzpp - klc
244             dpktCP(ikl,k) = PdTa_1(ikl,klc) /ExnrDY(ikl,k)
245             dqv_CP(ikl,k) = PdQa_1(ikl,klc)
246             dqw_CP(ikl,k) = PdQw_1(ikl,klc)
247             dqi_CP(ikl,k) = PdQi_1(ikl,klc)
248
249          END DO
250
251        END DO
252
253
254      END IF ! mod(iitCV0,jjtCV0).EQ.0   (UPDATE)                   THEN
255
256
257
258
259! Vertical Integrated Energy and Water Content
260! ============================================
261
262! #EW     irmx   =       i_x0
263! #EW     jrmx   =       j_y0
264! #EW     rr_max =       0.0
265
266        DO       ikl = 1,kcolp
267           i   = ii__AP(ikl)
268           j   = jj__AP(ikl)
269
270! #EW     enr0EW(ikl) = 0.0
271! #EW     wat0EW(ikl) = 0.0
272
273! #EW    DO k=1,mzp
274! #EW     temp_r      = pkt_DY(ikl,k)*ExnrDY(ikl,k)
275! #EW     enr0EW(ikl) = enr0EW(ikl)                                      &
276! #EW&                +(temp_r                                           &
277! #EW&                -(qw__CM(ikl,k)+qr__CM(ikl,k)) * Lv_CPd            &
278! #EW&                -(qi__CM(ikl,k)+qs__CM(ikl,k)) * Ls_CPd)*dsigmi(k)
279! #EW     wat0EW(ikl) = wat0EW(ikl)                                      &
280! #EW&                +(qv__DY(ikl,k)                                    &
281! #EW&                + qw__CM(ikl,k)+qr__CM(ikl,k)                      &
282! #EW&                + qi__CM(ikl,k)+qs__CM(ikl,k)          )*dsigmi(k)
283! #EW    END DO
284
285! #EW     enr0EW(ikl) = enr0EW(ikl) * psa_DY(  ikl)  * Grav_I
286! #EW     wat0EW(ikl) = wat0EW(ikl) * psa_DY(  ikl)  * Grav_I
287!  ..     wat0EW [m]    contains implicit factor 1.d3 [kPa-->Pa] /ro_Wat
288
289! #EW     energ0      = energ0 - enr0EW(ikl)
290! #EW     water0      = water0 - wat0EW(ikl)
291
292
293
294
295! Update of Mass Flux convective Tendencies
296! =========================================
297
298!       ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
299
300         IF                    (timeCP(ikl) .GT. 0.)                THEN
301
302
303!  Temporal tendencies on pkt_DY, qv__DY and rainCM
304!  ------------------------------------------------
305
306          DO k=1,mzp
307            dqv_CP(ikl,k) = min(dqv_CP(ikl,k),(qv__DY(ikl,k)-epsq)/dt__CP)
308            dqw_CP(ikl,k) = min(dqw_CP(ikl,k),(qw__CM(ikl,k)-epsq)/dt__CP)
309            dqi_CP(ikl,k) = min(dqi_CP(ikl,k),(qi__CM(ikl,k)-epsq)/dt__CP)
310          ENDDO
311
312            rainCM(ikl)   =     rainCM(ikl)  + drr_CP(ikl)        *dt__CP   
313            rainCP(ikl)   =     rainCP(ikl)  + drr_CP(ikl)        *dt__CP   
314
315            snowCM(ikl)   =     snowCM(ikl)  + dss_CP(ikl)        *dt__CP   
316            snowCP(ikl)   =     snowCP(ikl)  + dss_CP(ikl)        *dt__CP   
317
318         ELSE
319
320          DO k=1,mzp
321            dpktCP(ikl,k) =     0.00
322            dqv_CP(ikl,k) =     0.00
323            dqw_CP(ikl,k) =     0.00
324            dqi_CP(ikl,k) =     0.00
325          ENDDO
326
327
328         ENDIF          !      {timeCP(ikl) .gt. 0}
329
330            timeCP(ikl)   = max(timeCP(ikl) - dt__CP,zer0)
331!  ..       ^^^^ remaining time before the end of convection
332!       ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
333
334
335
336
337! Vertical Integrated Energy and Water Content
338! ============================================
339
340! #EW     enr1EW(ikl) = 0.0
341! #EW     wat1EW(ikl) = 0.0
342! #EW     watfEW(ikl) =-drr_CP(ikl)
343
344! #EW    DO k=1,mz
345! #EW     temp_r      = pkt_DY(ikl,k)*ExnrDY(ikl,k)
346! #EW     enr1EW(ikl) = enr1EW(ikl)                                      &
347! #EW&                +(temp_r                                           &
348! #EW&                -(qw__CM(ikl,k)+qr__CM(ikl,k)) *Lv_CPd             &
349! #EW&                -(qi__CM(ikl,k)+qs__CM(ikl,k)) *Ls_CPd)*dsigmi(k)
350! #EW     wat1EW(ikl) = wat1EW(ikl)                                      &
351! #EW&                +(qv__DY(ikl,k)                                    &
352! #EW&                + qw__CM(ikl,k)+qr__CM(ikl,k)                      &
353! #EW&                + qi__CM(ikl,k)+qs__CM(ikl,k)         )*dsigmi(k)
354! #EW    END DO
355
356! #EW     enr1EW(ikl) = enr1EW(ikl) * psa_DY(  ikl)  *Grav_I             &
357! #EW&                - drr_CP(ikl)                  *Lv_CPd
358! #EW     wat1EW(ikl) = wat1EW(ikl) * psa_DY(  ikl)  *Grav_I
359!  ..     wat1EW [m]    contains implicit factor 1.d3 [kPa-->Pa] /rhoWat
360
361! #EW     energ0      = energ0 + enr1EW(ikl)
362! #EW     water0      = water0 + wat1EW(ikl)
363! #EW     iter_0      = iter_0 + 1
364
365
366
367
368! Vertical Integrated Energy and Water Content: OUTPUT
369! ====================================================
370
371! #EW   IF (drr_CP(ikl).gt.rr_max)                                THEN
372! #EW       rr_max =       drr_CP(ikl)
373! #EW       irmx   =              i
374! #EW       jrmx   =                j
375! #EW   END IF
376
377        END DO
378
379! #EW   waterb = wat1EW(irmx,jrmx)-wat0EW(irmx,jrmx)-watfEW(irmx,jrmx)
380! #EW   write(6,606) it_EXP,enr0EW(irmx,jrmx),1.d3*wat0EW(irmx,jrmx),    &
381! #EW&            irmx,jrmx,enr1EW(irmx,jrmx),1.d3*wat1EW(irmx,jrmx),    &
382! #EW&                                        1.d3*watfEW(irmx,jrmx),    &
383! #EW&                                        1.d3*waterb           ,    &
384! #EW&                      energ0/iter_0    ,     water0/iter_0
385 606    format(i9,'  Before CVAj:  E0 =',f12.6,'  W0 = ',f9.6,           &
386     &    /,i5,i4,'  After  CVAj:  E1 =',f12.6,'  W1 = ',f9.6,           &
387     &                                         '  W Flux =',f9.6,        &
388     &                                         '  Div(W) =',e9.3,        &
389     &       /,9x,'         Mean   dE =',f12.9,'  dW = ',e9.3)
390
391
392
393
394! Incrementation Step Nb
395! ======================
396
397      iitCV0 = iitCV0 + 1
398
399
400!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
401!                                                                                !
402!  DE-ALLOCATION
403!  =============
404
405      IF (FlagDALLOC)                                      THEN !
406          deallocate  ( wa_ANA )                                                 !  ANAbatic       Wind        Speed     [m/s]
407      END IF
408
409! Bug corrigé par Martin:
410!      IF (it_RUN.EQ.1 .OR. FlagDALLOC)                                      THEN !
411!          deallocate  ( wa_ANA )                                                 !  ANAbatic       Wind        Speed     [m/s]
412!      END IF
413!                                                                                !
414!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
415
416
417
418      return
419      end subroutine PHY_Atm_CP_RUN
Note: See TracBrowser for help on using the repository browser.