source: LMDZ5/trunk/libf/phymar/PHY_Atm_AT_RUN.f90 @ 5412

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

Inclusion de la physique de MAR


Integration of MAR physics

File size: 13.8 KB
Line 
1      subroutine PHY_Atm_AT_RUN(FlagSV_KzT,FlagCM)
2
3!------------------------------------------------------------------------------+
4!                                                         Sat 22-Jun-2013  MAR |
5!   MAR          PHY_Atm_AT_RUN                                                |
6!     subroutine PHY_Atm_AT_RUN performs   Turbulent Vertical Diffusion Scheme |
7!                                                                              |
8!     version 3.p.4.1 created by H. Gallee,               Wed 13-Mar-2013      |
9!           Last Modification by H. Gallee,               Sat 22-Jun-2013      |
10!                                                                              |
11!------------------------------------------------------------------------------+
12
13      use Mod_Real
14      use Mod_PHY____dat
15      use Mod_PHY____grd
16      use Mod_PHY_AT_ctr
17      use Mod_PHY_AT_grd
18      use Mod_PHY_AT_kkl
19      use Mod_SISVAT_gpt
20
21
22      IMPLICIT NONE
23
24
25      logical    ::    FlagSV_KzT
26      logical    ::    FlagCM
27
28      integer    ::    i,     j,     k,     ikl
29
30
31
32
33! Assignation    of Mod_PHY_AT_kkl
34! ================================
35
36        DO ikl = 1,kcolp
37
38          i = ii__AP(ikl)
39          j = jj__AP(ikl)
40
41          LMO_AT(ikl)   = LMO_SV_gpt(ikl)
42
43        ENDDO
44
45
46
47
48! Turbulent Kinetic Energy
49! ========================
50
51      IF (AT_TKE)                                                   THEN
52
53!              **************
54        CALL   PHY_vdfTKE_RUN
55        CALL   PHY_genTKE_RUN
56!              **************
57
58      ELSE
59         DO ikl = 1,kcolp
60         DO   k = 1,mzp
61           TKE_AT(ikl,k)=    eps6
62           eps_AT(ikl,k)=    epsn
63         ENDDO
64         ENDDO
65      END IF
66
67
68
69
70! Turbulent Vertical Diffusion
71! ============================
72 
73          NewAAT =  .TRUE.
74          SchmAT =  'u'   
75!                    ******
76          CALL       Atm_AT
77!                    ******
78
79
80          NewAAT =  .FALSE.
81          SchmAT =  'v'   
82!                    ******
83          CALL       Atm_AT
84!                    ******
85
86
87      IF (.NOT.FlagSV_KzT)                                          THEN
88
89          SchmAT =  'T'   
90!                    ******
91          CALL       Atm_AT
92!                    ******
93      ELSE
94
95          DO k  =1,mzp
96          DO ikl=1,kcolp
97            dpktAT(ikl,k) = 0.0000
98          ENDDO
99          ENDDO
100
101      END IF
102
103
104          SchmAT =  'q'   
105!                    ******
106          CALL       Atm_AT
107!                    ******
108
109
110      IF (FlagCM)                                                   THEN
111
112          SchmAT =  'w'   
113!                    ******
114          CALL       Atm_AT
115!                    ******
116
117
118          SchmAT =  'i'   
119!                    ******
120          CALL       Atm_AT
121!                    ******
122
123
124          SchmAT =  'C'   
125!                    ******
126          CALL       Atm_AT
127!                    ******
128
129
130          SchmAT =  's'   
131!                    ******
132          CALL       Atm_AT
133!                    ******
134
135
136          SchmAT =  'r'   
137!                    ******
138          CALL       Atm_AT
139!                    ******
140
141
142      END IF
143
144
145      end subroutine PHY_Atm_AT_RUN
146
147
148
149
150      subroutine Atm_AT
151
152!------------------------------------------------------------------------------+
153!                                                         Sat 22-Jun-2013  MAR |
154!   MAR          Atm_AT                                                        |
155!     subroutine Atm_AT solves Turbulent Vertical Diffusion Scheme             |
156!                                                                              |
157!     version 3.p.4.1 created by H. Gallee,               Wed 13-Mar-2013      |
158!           Last Modification by H. Gallee,               Sat 22-Jun-2013      |
159!                                                                              |
160!------------------------------------------------------------------------------+
161
162      use Mod_Real
163      use Mod_PHY____dat
164      use Mod_PHY____grd
165      use Mod_PHY_AT_grd
166      use Mod_PHY_AT_kkl
167      use Mod_PHY_CM_kkl
168      use Mod_PHY_DY_kkl
169      use Mod_SISVAT_gpt
170
171
172
173      IMPLICIT NONE
174
175
176
177      integer           ::       i,     j,     k,     ikl
178      real(kind=real8)  ::  psa_sq
179
180
181
182
183!     ++++++++++++++++
184      DO ikl = 1,kcolp
185!     ++++++++++++++++
186
187        i = ii__AP(ikl)
188        j = jj__AP(ikl)
189
190
191
192! Tri-Diagonal Matrix: Variable
193! =============================
194
195!  Wind Speed
196!  ----------
197
198        IF      (schmAT.EQ.'a')                                     THEN
199          DO k=1,mzp
200            Kz__AT(ikl,k) = 1.0
201            var_AT(ikl,k) = ua__DY(ikl,k)
202          ENDDO
203             k=  mzpp
204            var_AT(ikl,k) = 0.
205
206        ELSE IF (schmAT.EQ.'b')                                     THEN
207          DO k=1,mzp
208            Kz__AT(ikl,k) = 1.0
209            var_AT(ikl,k) = va__DY(ikl,k)
210          ENDDO
211             k=  mzpp
212            var_AT(ikl,k) = 0.
213
214        ELSE IF (schmAT.EQ.'u')                                     THEN
215          DO k=1,mzp
216            Kz__AT(ikl,k) = Kzm_AT(ikl,k)
217            var_AT(ikl,k) = ua__DY(ikl,k)
218          ENDDO
219             k=  mzpp
220            var_AT(ikl,k) = 0.
221
222        ELSE IF (schmAT.EQ.'v')                                     THEN
223          DO k=1,mzp
224            Kz__AT(ikl,k) = Kzm_AT(ikl,k)
225            var_AT(ikl,k) = va__DY(ikl,k)
226          ENDDO
227             k=  mzpp
228            var_AT(ikl,k) = 0.
229
230
231!  Temperature
232!  -----------
233
234        ELSE IF (schmAT.EQ.'T')                                     THEN
235          DO k=1,mzp
236            Kz__AT(ikl,k) = Kzh_AT(ikl,k)
237            var_AT(ikl,k) = pkt_DY(ikl,k)
238          ENDDO
239             k=  mzpp
240            var_AT(ikl,k) = pkt_DY(ikl,k)
241
242
243!  Water Species
244!  -------------
245
246        ELSE IF (schmAT.EQ.'q')                                     THEN
247          DO k=1,mzp
248            Kz__AT(ikl,k) = Kzh_AT(ikl,k)
249            var_AT(ikl,k) = qv__DY(ikl,k)
250          ENDDO
251             k=  mzpp
252            var_AT(ikl,k) = qv__DY(ikl,mzp)                            &
253     &                    -(Z___DY(ikl,mzp)-Z___DY(ikl,mzpp))          &
254     &                    * uqs_SV_gpt(ikl)    /Kzh_AT(ikl,mzp )
255            var_AT(ikl,k) =            max(epsq,var_AT(ikl,k   ))
256            var_AT(ikl,k) =            min(     var_AT(ikl,k   ),qvsiCM(ikl,k))
257            qv__DY(ikl,mzpp) =                  var_AT(ikl,k   )
258
259        ELSE IF (schmAT.EQ.'w')                                     THEN
260          DO k=1,mzp
261            Kz__AT(ikl,k) = Kzh_AT(ikl,k)
262            var_AT(ikl,k) = qw__CM(ikl,k)
263          ENDDO
264             k=  mzpp
265            var_AT(ikl,k) = qw__CM(ikl,mzp)
266
267        ELSE IF (schmAT.EQ.'i')                                     THEN
268          DO k=1,mzp
269            Kz__AT(ikl,k) = Kzh_AT(ikl,k)
270            var_AT(ikl,k) = qi__CM(ikl,k)
271          ENDDO
272             k=  mzpp
273            var_AT(ikl,k) = qi__CM(ikl,mzp)
274
275        ELSE IF (schmAT.EQ.'C')                                     THEN
276          DO k=1,mzp
277            Kz__AT(ikl,k) = Kzh_AT(ikl,k)
278            var_AT(ikl,k) = CCNiCM(ikl,k)
279          ENDDO
280             k=  mzpp
281            var_AT(ikl,k) = CCNiCM(ikl,mzp)
282
283        ELSE IF (schmAT.EQ.'s')                                     THEN
284          DO k=1,mzp
285            Kz__AT(ikl,k) = Kzh_AT(ikl,k)
286            var_AT(ikl,k) = qs__CM(ikl,k)
287          ENDDO
288             k=  mzpp
289            var_AT(ikl,k) = qs__CM(ikl,k)
290
291        ELSE IF (schmAT.EQ.'r')                                     THEN
292          DO k=1,mzp
293            Kz__AT(ikl,k) = Kzh_AT(ikl,k)
294            var_AT(ikl,k) = qr__CM(ikl,k)
295          ENDDO
296             k=  mzpp
297            var_AT(ikl,k) = qr__CM(ikl,mzp)
298
299        ELSE
300            write(6,*)
301            write(6,*) ' schmAT = ',schmAT
302            write(6,*) ' Inadequate..., NOR for Ekman Spiral, NOR for momentum, NOR for scalars'
303            write(6,*) '                                                    '
304
305!                **************
306            call MAR________BUG
307!                **************
308            STOP
309        END IF
310
311
312
313! Tri-Diagonal Matrix Coefficients
314! ================================
315
316        IF (NewAAT)                                                 THEN
317            psa_sq          = psa_DY(ikl)     *psa_DY(ikl)   
318            Ac__AT(ikl,mzp) = Ac0_AT(mzp)                                   &
319     &                      * roa_DY(ikl,mzp) *roamDY(ikl,mzp)              &
320     &                      / psa_sq
321          DO k=mzp-1,1,-1
322            Ac__AT(ikl,k  ) = Ac0_AT(k  )                                   &
323     &                      * roa_DY(ikl,k  ) *roamDY(ikl,k  )              &
324     &                      / psa_sq
325            Cc__AT(ikl,k+1) = Ac__AT(ikl,k)                                 &
326     &                      * roa_DY(ikl,k+1) /roa_DY(ikl,k  )
327          ENDDO
328        ENDIF
329
330            A___AT(ikl,mzp) =           Ac__AT(ikl,mzp) * Kz__AT(ikl,mzp)
331          DO k=mzp-1,1,-1
332            A___AT(ikl,k  ) =           Ac__AT(ikl,k  ) * Kz__AT(ikl,k  )
333            C___AT(ikl,k+1) =           Cc__AT(ikl,k+1) * Kz__AT(ikl,k  )
334          ENDDO
335          DO k=    2,mzp
336            B___AT(ikl,k  ) =  1.0000 - A___AT(ikl,k  ) - C___AT(ikl,k  )
337          ENDDO
338            B___AT(ikl,1  ) =  1.0000 - A___AT(ikl,1  )
339            C___AT(ikl,1  ) =  0.0000
340
341            D___AT(ikl,1  ) =                                               &
342     &      A___AT(ikl,1  ) *  a_b_AT *(var_AT(ikl,1  ) - var_AT(ikl,2  ))
343
344             k=      mzp
345            D___AT(ikl,k  ) =                                               &
346     &      A___AT(ikl,k  ) * (a_b_AT * var_AT(ikl,k  ) - var_AT(ikl,k+1)   &
347     &                                                  / betaAT)           &
348     &    - C___AT(ikl,k  ) *  a_b_AT *(var_AT(ikl,k-1) - var_AT(ikl,k  ))  &
349     &    +                                               var_AT(ikl,k  )
350            A___AT(ikl,k  ) =  0.0000
351          DO k=    2,mzp-1
352            D___AT(ikl,k  ) =                                               &
353     &      A___AT(ikl,k  ) *  a_b_AT *(var_AT(ikl,k  ) - var_AT(ikl,k+1))  &
354     &    - C___AT(ikl,k  ) *  a_b_AT *(var_AT(ikl,k-1) - var_AT(ikl,k  ))  &
355     &    +                                               var_AT(ikl,k  )
356          ENDDO
357            D___AT(ikl,1  ) =                                               &
358     &      A___AT(ikl,1  ) *  a_b_AT *(var_AT(ikl,1  ) - var_AT(ikl,2  ))  &
359     &    +                                               var_AT(ikl,1  )
360
361
362
363! Tri-Diagonal Matrix Resolution (Gaussian Elimination / Thomas Algorithm)
364! ==============================
365
366! Forward  Sweep
367! --------------
368
369          P___AT(1) =  B___AT(ikl,1)
370          Q___AT(1) = -A___AT(ikl,1) / P___AT(1)
371        DO k=2,mzp
372          P___AT(k) =  C___AT(ikl,k) * Q___AT(k-1)   + B___AT(ikl,k)
373          Q___AT(k) = -A___AT(ikl,k) / P___AT(k)
374        ENDDO
375          X___AT(1) =  D___AT(ikl,1) / P___AT(1)
376        DO k=2,mzp
377          X___AT(k) = (D___AT(ikl,k) - C___AT(ikl,k) * X___AT(k-1)) / P___AT(k)
378        END DO
379
380
381! Backward Sweep
382! --------------
383
384        DO k=  mzp-1,1,-1
385          X___AT(k) =  Q___AT(k)     * X___AT(k+1)   + X___AT(k)
386        END DO
387
388
389
390! Elimination of too small values
391! -------------------------------
392
393        IF (SchmAT.EQ.'w'.OR.SchmAT.EQ.'i'.OR.SchmAT.EQ.'C'.OR.        &
394     &      SchmAT.EQ.'s'.OR.SchmAT.EQ.'r')                         THEN
395          DO k=  mzp  ,1,-1
396            IF (X___AT(k).NE.0.0000 .AND. X___AT(k).LT.1.e-18)         &
397     &          X___AT(k) =  0.0000
398          END DO
399        END IF
400
401
402
403
404! Tendency
405! ========
406
407!  Wind Speed
408!  ----------
409
410        IF      (schmAT.EQ.'a')                                     THEN
411          DO k=1,mzp
412            dua_AT(ikl,k) =(X___AT(k) - ua__DY(ikl,k))  / dt__AT
413          ENDDO
414        ELSE IF (schmAT.EQ.'b')                                     THEN
415          DO k=1,mzp
416            dva_AT(ikl,k) =(X___AT(k) - va__DY(ikl,k))  / dt__AT
417          ENDDO
418        ELSE IF (schmAT.EQ.'u')                                     THEN
419          DO k=1,mzp
420            dua_AT(ikl,k) =(X___AT(k) - ua__DY(ikl,k))  / dt__AT
421          ENDDO
422        ELSE IF (schmAT.EQ.'v')                                     THEN
423          DO k=1,mzp
424            dva_AT(ikl,k) =(X___AT(k) - va__DY(ikl,k))  / dt__AT
425          ENDDO
426
427
428!  Temperature
429!  -----------
430
431        ELSE IF (schmAT.EQ.'T')                                     THEN
432          DO k=1,mzp
433            dpktAT(ikl,k) =(X___AT(k) - pkt_DY(ikl,k))  / dt__AT
434          ENDDO
435
436
437!  Water Species
438!  -------------
439
440        ELSE IF (schmAT.EQ.'q')                                     THEN
441          DO k=1,mzp
442            dqv_AT(ikl,k) =(X___AT(k) - qv__DY(ikl,k))  / dt__AT
443          ENDDO
444
445        ELSE IF (schmAT.EQ.'w')                                     THEN
446          DO k=1,mzp
447            dqw_AT(ikl,k) =(X___AT(k) - qw__CM(ikl,k))  / dt__AT
448          ENDDO
449
450        ELSE IF (schmAT.EQ.'i')                                     THEN
451          DO k=1,mzp
452            dqi_AT(ikl,k) =(X___AT(k) - qi__CM(ikl,k))  / dt__AT
453          ENDDO
454
455        ELSE IF (schmAT.EQ.'C')                                     THEN
456          DO k=1,mzp
457            dCi_AT(ikl,k) =(X___AT(k) - CCNiCM(ikl,k))  / dt__AT
458          ENDDO
459
460        ELSE IF (schmAT.EQ.'s')                                     THEN
461          DO k=1,mzp
462            dqs_AT(ikl,k) =(X___AT(k) - qs__CM(ikl,k))  / dt__AT
463          ENDDO
464
465        ELSE IF (schmAT.EQ.'r')                                     THEN
466          DO k=1,mzp
467            dqr_AT(ikl,k) =(X___AT(k) - qr__CM(ikl,k))  / dt__AT
468          ENDDO
469
470        ELSE
471            write(6,*)
472            write(6,*) ' schmAT = ',schmAT
473            write(6,*) ' Inadequate..., NOR for Ekman Spiral, NOR for momentum, NOR for scalars'
474            write(6,*) '                                                    '
475
476!                **************
477            call MAR________BUG
478!                **************
479            STOP
480        END IF
481
482
483!     ++++++
484      END DO
485!     ++++++
486
487
488
489
490      end subroutine Atm_AT
Note: See TracBrowser for help on using the repository browser.