source: LMDZ6/trunk/libf/phylmd/climb_wind_mod.f90 @ 5279

Last change on this file since 5279 was 5274, checked in by abarral, 9 months ago

Replace yomcst.h by existing module
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 15.4 KB
Line 
1!
2MODULE climb_wind_mod
3!
4! Module to solve the verctical diffusion of the wind components "u" and "v".
5!
6  USE dimphy
7
8  IMPLICIT NONE
9
10  SAVE
11  PRIVATE
12 
13  REAL, DIMENSION(:),   ALLOCATABLE  :: alf1, alf2
14  !$OMP THREADPRIVATE(alf1,alf2)
15  REAL, DIMENSION(:,:), ALLOCATABLE  :: Kcoefm
16  !$OMP THREADPRIVATE(Kcoefm)
17  REAL, DIMENSION(:,:), ALLOCATABLE  :: Ccoef_U, Dcoef_U
18  !$OMP THREADPRIVATE(Ccoef_U, Dcoef_U)
19  REAL, DIMENSION(:,:), ALLOCATABLE  :: Ccoef_V, Dcoef_V
20  !$OMP THREADPRIVATE(Ccoef_V, Dcoef_V)
21  REAL, DIMENSION(:), ALLOCATABLE   :: Acoef_U, Bcoef_U
22  !$OMP THREADPRIVATE(Acoef_U, Bcoef_U)
23  REAL, DIMENSION(:), ALLOCATABLE   :: Acoef_V, Bcoef_V
24  !$OMP THREADPRIVATE(Acoef_V, Bcoef_V)
25  LOGICAL                            :: firstcall=.TRUE.
26  !$OMP THREADPRIVATE(firstcall)
27
28 
29  PUBLIC :: climb_wind_down, climb_wind_up
30
31CONTAINS
32!
33!****************************************************************************************
34!
35  SUBROUTINE climb_wind_init
36
37    INTEGER             :: ierr
38    CHARACTER(len = 20) :: modname = 'climb_wind_init'   
39
40!****************************************************************************************
41! Allocation of global module variables
42!
43!****************************************************************************************
44
45    ALLOCATE(alf1(klon), stat=ierr)
46    IF (ierr /= 0) CALL abort_physic(modname,'Pb in allocate alf1',1)
47
48    ALLOCATE(alf2(klon), stat=ierr)
49    IF (ierr /= 0) CALL abort_physic(modname,'Pb in allocate alf2',1)
50
51    ALLOCATE(Kcoefm(klon,klev), stat=ierr)
52    IF (ierr /= 0) CALL abort_physic(modname,'Pb in allocate Kcoefm',1)
53
54    ALLOCATE(Ccoef_U(klon,klev), stat=ierr)
55    IF (ierr /= 0) CALL abort_physic(modname,'Pb in allocate Ccoef_U',1)
56
57    ALLOCATE(Dcoef_U(klon,klev), stat=ierr)
58    IF (ierr /= 0) CALL abort_physic(modname,'Pb in allocation Dcoef_U',1)
59
60    ALLOCATE(Ccoef_V(klon,klev), stat=ierr)
61    IF (ierr /= 0) CALL abort_physic(modname,'Pb in allocation Ccoef_V',1)
62
63    ALLOCATE(Dcoef_V(klon,klev), stat=ierr)
64    IF (ierr /= 0) CALL abort_physic(modname,'Pb in allocation Dcoef_V',1)
65
66    ALLOCATE(Acoef_U(klon), Bcoef_U(klon), Acoef_V(klon), Bcoef_V(klon), STAT=ierr)
67    IF ( ierr /= 0 )  PRINT*,' pb in allloc Acoef_U and Bcoef_U, ierr=', ierr
68
69    firstcall=.FALSE.
70
71  END SUBROUTINE climb_wind_init
72!
73!****************************************************************************************
74!
75  SUBROUTINE climb_wind_down(knon, dtime, coef_in, pplay, paprs, temp, delp, u_old, v_old, &
76!!! nrlmd le 02/05/2011
77       Ccoef_U_out, Ccoef_V_out, Dcoef_U_out, Dcoef_V_out, &
78       Kcoef_m_out, alf_1_out, alf_2_out, &
79!!!
80       Acoef_U_out, Acoef_V_out, Bcoef_U_out, Bcoef_V_out)
81!
82! This routine calculates for the wind components u and v,
83! recursivly the coefficients C and D in equation
84! X(k) = C(k) + D(k)*X(k-1), X=[u,v], k=[1,klev] is the vertical layer.
85!
86!
87USE yomcst_mod_h, ONLY: RPI, RCLUM, RHPLA, RKBOL, RNAVO                   &
88          , RDAY, REA, REPSM, RSIYEA, RSIDAY, ROMEGA                  &
89          , R_ecc, R_peri, R_incl                                      &
90          , RA, RG, R1SA                                         &
91          , RSIGMA                                                     &
92          , R, RMD, RMV, RD, RV, RCPD                    &
93          , RMO3, RMCO2, RMC, RMCH4, RMN2O, RMCFC11, RMCFC12        &
94          , RCPV, RCVD, RCVV, RKAPPA, RETV, eps_w                    &
95          , RCW, RCS                                                 &
96          , RLVTT, RLSTT, RLMLT, RTT, RATM                           &
97          , RESTT, RALPW, RBETW, RGAMW, RALPS, RBETS, RGAMS            &
98          , RALPD, RBETD, RGAMD
99! Input arguments
100!****************************************************************************************
101    INTEGER, INTENT(IN)                      :: knon
102    REAL, INTENT(IN)                         :: dtime
103    REAL, DIMENSION(klon,klev), INTENT(IN)   :: coef_in
104    REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay ! pres au milieu de couche (Pa)
105    REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! pression a inter-couche (Pa)
106    REAL, DIMENSION(klon,klev), INTENT(IN)   :: temp  ! temperature
107    REAL, DIMENSION(klon,klev), INTENT(IN)   :: delp
108    REAL, DIMENSION(klon,klev), INTENT(IN)   :: u_old
109    REAL, DIMENSION(klon,klev), INTENT(IN)   :: v_old
110
111! Output arguments
112!****************************************************************************************
113    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef_U_out
114    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef_V_out
115    REAL, DIMENSION(klon), INTENT(OUT)       :: Bcoef_U_out
116    REAL, DIMENSION(klon), INTENT(OUT)       :: Bcoef_V_out
117
118!!! nrlmd le 02/05/2011
119    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Ccoef_U_out
120    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Ccoef_V_out
121    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Dcoef_U_out
122    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Dcoef_V_out
123    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Kcoef_m_out
124    REAL, DIMENSION(klon), INTENT(OUT)       :: alf_1_out
125    REAL, DIMENSION(klon), INTENT(OUT)       :: alf_2_out
126!!!
127
128! Local variables
129!****************************************************************************************
130    REAL, DIMENSION(klon)                    :: u1lay, v1lay
131    INTEGER                                  :: k, i
132
133! Include
134!****************************************************************************************
135    INCLUDE "compbl.h"
136
137!****************************************************************************************
138! Initialize module
139    IF (firstcall) CALL climb_wind_init
140
141!****************************************************************************************
142! Calculate the coefficients C and D in : u(k) = C(k) + D(k)*u(k-1)
143!
144!****************************************************************************************
145! - Define alpha (alf1 and alf2)
146    alf1(:) = 1.0
147    alf2(:) = 1.0 - alf1(:)
148
149! - Calculate the coefficients K
150    Kcoefm(:,:) = 0.0
151    DO k = 2, klev
152       DO i=1,knon
153          Kcoefm(i,k) = coef_in(i,k)*RG*RG*dtime/(pplay(i,k-1)-pplay(i,k)) &
154               *(paprs(i,k)*2/(temp(i,k)+temp(i,k-1))/RD)**2
155       END DO
156    END DO
157
158! - Calculate the coefficients C and D, component "u"
159    CALL calc_coef(knon, Kcoefm(:,:), delp(:,:), &
160         u_old(:,:), alf1(:), alf2(:),  &
161         Ccoef_U(:,:), Dcoef_U(:,:), Acoef_U(:), Bcoef_U(:))
162
163! - Calculate the coefficients C and D, component "v"
164    CALL calc_coef(knon, Kcoefm(:,:), delp(:,:), &
165         v_old(:,:), alf1(:), alf2(:),  &
166         Ccoef_V(:,:), Dcoef_V(:,:), Acoef_V(:), Bcoef_V(:))
167
168!****************************************************************************************
169! 6)
170! Return the first layer in output variables
171!
172!****************************************************************************************
173    Acoef_U_out = Acoef_U
174    Bcoef_U_out = Bcoef_U
175    Acoef_V_out = Acoef_V
176    Bcoef_V_out = Bcoef_V
177
178!****************************************************************************************
179! 7)
180! If Pbl is split, return also the other layers in output variables
181!
182!****************************************************************************************
183!!! jyg le 07/02/2012
184!!jyg       IF (mod(iflag_pbl_split,2) .eq.1) THEN
185       IF (mod(iflag_pbl_split,10) .ge.1) THEN
186!!! nrlmd le 02/05/2011
187    DO k= 1, klev
188      DO i= 1, klon
189        Ccoef_U_out(i,k) = Ccoef_U(i,k)
190        Ccoef_V_out(i,k) = Ccoef_V(i,k)
191        Dcoef_U_out(i,k) = Dcoef_U(i,k)
192        Dcoef_V_out(i,k) = Dcoef_V(i,k)
193        Kcoef_m_out(i,k) = Kcoefm(i,k)
194      ENDDO
195    ENDDO
196    DO i= 1, klon
197      alf_1_out(i)   = alf1(i)
198      alf_2_out(i)   = alf2(i)
199    ENDDO
200!!!     
201       ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
202!!!
203
204  END SUBROUTINE climb_wind_down
205!
206!****************************************************************************************
207!
208  SUBROUTINE calc_coef(knon, Kcoef, delp, X, alfa1, alfa2, Ccoef, Dcoef, Acoef, Bcoef)
209!
210! Find the coefficients C and D in fonction of alfa, K and delp
211USE yomcst_mod_h, ONLY: RPI, RCLUM, RHPLA, RKBOL, RNAVO                   &
212          , RDAY, REA, REPSM, RSIYEA, RSIDAY, ROMEGA                  &
213          , R_ecc, R_peri, R_incl                                      &
214          , RA, RG, R1SA                                         &
215          , RSIGMA                                                     &
216          , R, RMD, RMV, RD, RV, RCPD                    &
217          , RMO3, RMCO2, RMC, RMCH4, RMN2O, RMCFC11, RMCFC12        &
218          , RCPV, RCVD, RCVV, RKAPPA, RETV, eps_w                    &
219          , RCW, RCS                                                 &
220          , RLVTT, RLSTT, RLMLT, RTT, RATM                           &
221          , RESTT, RALPW, RBETW, RGAMW, RALPS, RBETS, RGAMS            &
222          , RALPD, RBETD, RGAMD
223! Input arguments
224!****************************************************************************************
225    INTEGER, INTENT(IN)                      :: knon
226    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Kcoef, delp
227    REAL, DIMENSION(klon,klev), INTENT(IN)   :: X
228    REAL, DIMENSION(klon), INTENT(IN)        :: alfa1, alfa2
229
230! Output arguments
231!****************************************************************************************
232    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef, Bcoef
233    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Ccoef, Dcoef
234 
235! local variables
236!****************************************************************************************
237    INTEGER                                  :: k, i
238    REAL                                     :: buf
239    !****************************************************************************************
240!
241
242! Calculate coefficients C and D at top level, k=klev
243!
244    Ccoef(:,:) = 0.0
245    Dcoef(:,:) = 0.0
246
247    DO i = 1, knon
248       buf = delp(i,klev) + Kcoef(i,klev)
249
250       Ccoef(i,klev) = X(i,klev)*delp(i,klev)/buf
251       Dcoef(i,klev) = Kcoef(i,klev)/buf
252    END DO
253   
254!
255! Calculate coefficients C and D at top level (klev-1) <= k <= 2
256!
257    DO k=(klev-1),2,-1
258       DO i = 1, knon
259          buf = delp(i,k) + Kcoef(i,k) + Kcoef(i,k+1)*(1.-Dcoef(i,k+1))
260         
261          Ccoef(i,k) = (X(i,k)*delp(i,k) + Kcoef(i,k+1)*Ccoef(i,k+1))/buf
262          Dcoef(i,k) = Kcoef(i,k)/buf
263       END DO
264    END DO
265
266!
267! Calculate coeffiecent A and B at surface
268!
269    DO i = 1, knon
270       buf = delp(i,1) + Kcoef(i,2)*(1-Dcoef(i,2))
271       Acoef(i) = (X(i,1)*delp(i,1) + Kcoef(i,2)*Ccoef(i,2))/buf
272       Bcoef(i) = -RG/buf
273    END DO
274
275  END SUBROUTINE calc_coef
276!
277!****************************************************************************************
278!
279
280  SUBROUTINE climb_wind_up(knon, dtime, u_old, v_old, flx_u1, flx_v1,  &
281!!! nrlmd le 02/05/2011
282       Acoef_U_in, Acoef_V_in, Bcoef_U_in, Bcoef_V_in, &
283       Ccoef_U_in, Ccoef_V_in, Dcoef_U_in, Dcoef_V_in, &
284       Kcoef_m_in, &
285!!!
286       flx_u_new, flx_v_new, d_u_new, d_v_new)
287!
288! Diffuse the wind components from the surface layer and up to the top layer.
289! Coefficents A, B, C and D are known from before. Start values for the diffusion are the
290! momentum fluxes at surface.
291!
292! u(k=1) = A + B*flx*dtime
293! u(k)   = C(k) + D(k)*u(k-1)  [2 <= k <= klev]
294!
295!****************************************************************************************
296USE yomcst_mod_h, ONLY: RPI, RCLUM, RHPLA, RKBOL, RNAVO                   &
297          , RDAY, REA, REPSM, RSIYEA, RSIDAY, ROMEGA                  &
298          , R_ecc, R_peri, R_incl                                      &
299          , RA, RG, R1SA                                         &
300          , RSIGMA                                                     &
301          , R, RMD, RMV, RD, RV, RCPD                    &
302          , RMO3, RMCO2, RMC, RMCH4, RMN2O, RMCFC11, RMCFC12        &
303          , RCPV, RCVD, RCVV, RKAPPA, RETV, eps_w                    &
304          , RCW, RCS                                                 &
305          , RLVTT, RLSTT, RLMLT, RTT, RATM                           &
306          , RESTT, RALPW, RBETW, RGAMW, RALPS, RBETS, RGAMS            &
307          , RALPD, RBETD, RGAMD
308! Input arguments
309!****************************************************************************************
310    INTEGER, INTENT(IN)                     :: knon
311    REAL, INTENT(IN)                        :: dtime
312    REAL, DIMENSION(klon,klev), INTENT(IN)  :: u_old
313    REAL, DIMENSION(klon,klev), INTENT(IN)  :: v_old
314    REAL, DIMENSION(klon), INTENT(IN)       :: flx_u1, flx_v1 ! momentum flux
315
316!!! nrlmd le 02/05/2011
317    REAL, DIMENSION(klon), INTENT(IN)       :: Acoef_U_in,Acoef_V_in, Bcoef_U_in, Bcoef_V_in
318    REAL, DIMENSION(klon,klev), INTENT(IN)  :: Ccoef_U_in, Ccoef_V_in, Dcoef_U_in, Dcoef_V_in
319    REAL, DIMENSION(klon,klev), INTENT(IN)  :: Kcoef_m_in
320!!!
321
322! Output arguments
323!****************************************************************************************
324    REAL, DIMENSION(klon,klev), INTENT(OUT) :: flx_u_new, flx_v_new
325    REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_u_new, d_v_new
326
327! Local variables
328!****************************************************************************************
329    REAL, DIMENSION(klon,klev)              :: u_new, v_new
330    INTEGER                                 :: k, i
331
332! Include
333!****************************************************************************************
334    INCLUDE "compbl.h"
335   
336!
337!****************************************************************************************
338
339!!! jyg le 07/02/2012
340!!jyg       IF (mod(iflag_pbl_split,2) .eq.1) THEN
341       IF (mod(iflag_pbl_split,10) .ge.1) THEN
342!!! nrlmd le 02/05/2011
343    DO i = 1, knon
344      Acoef_U(i)=Acoef_U_in(i)
345      Acoef_V(i)=Acoef_V_in(i)
346      Bcoef_U(i)=Bcoef_U_in(i)
347      Bcoef_V(i)=Bcoef_V_in(i)
348    ENDDO
349    DO k = 1, klev
350      DO i = 1, knon
351        Ccoef_U(i,k)=Ccoef_U_in(i,k)
352        Ccoef_V(i,k)=Ccoef_V_in(i,k)
353        Dcoef_U(i,k)=Dcoef_U_in(i,k)
354        Dcoef_V(i,k)=Dcoef_V_in(i,k)
355        Kcoefm(i,k)=Kcoef_m_in(i,k)
356      ENDDO
357    ENDDO
358!!!
359       ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
360!!!
361
362! Niveau 1
363    DO i = 1, knon
364       u_new(i,1) = Acoef_U(i) + Bcoef_U(i)*flx_u1(i)*dtime
365       v_new(i,1) = Acoef_V(i) + Bcoef_V(i)*flx_v1(i)*dtime
366    END DO
367
368! Niveau 2 jusqu'au sommet klev
369    DO k = 2, klev
370       DO i=1, knon
371          u_new(i,k) = Ccoef_U(i,k) + Dcoef_U(i,k) * u_new(i,k-1)
372          v_new(i,k) = Ccoef_V(i,k) + Dcoef_V(i,k) * v_new(i,k-1)
373       END DO
374    END DO
375
376!****************************************************************************************
377! Calcul flux
378!
379!== flux_u/v est le flux de moment angulaire (positif vers bas)
380!== dont l'unite est: (kg m/s)/(m**2 s)
381!
382!****************************************************************************************
383!
384    flx_u_new(:,:) = 0.0
385    flx_v_new(:,:) = 0.0
386
387    flx_u_new(1:knon,1)=flx_u1(1:knon)
388    flx_v_new(1:knon,1)=flx_v1(1:knon)
389
390! Niveau 2->klev
391    DO k = 2, klev
392       DO i = 1, knon
393          flx_u_new(i,k) = Kcoefm(i,k)/RG/dtime * &
394               (u_new(i,k)-u_new(i,k-1))
395         
396          flx_v_new(i,k) = Kcoefm(i,k)/RG/dtime * &
397               (v_new(i,k)-v_new(i,k-1))
398       END DO
399    END DO
400
401!****************************************************************************************
402! Calcul tendances
403!
404!****************************************************************************************
405    d_u_new(:,:) = 0.0
406    d_v_new(:,:) = 0.0
407    DO k = 1, klev
408       DO i = 1, knon
409          d_u_new(i,k) = u_new(i,k) - u_old(i,k)
410          d_v_new(i,k) = v_new(i,k) - v_old(i,k)
411       END DO
412    END DO
413
414  END SUBROUTINE climb_wind_up
415!
416!****************************************************************************************
417!
418END MODULE climb_wind_mod
Note: See TracBrowser for help on using the repository browser.