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

Last change on this file since 5300 was 5296, checked in by abarral, 4 weeks ago

Turn compbl.h into a module
Move calcul_REGDYN.h to obsolete
Create phys_constants_mod.f90
  • 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: 12.9 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
88USE compbl_mod_h
89! Input arguments
90!****************************************************************************************
91    INTEGER, INTENT(IN)                      :: knon
92    REAL, INTENT(IN)                         :: dtime
93    REAL, DIMENSION(klon,klev), INTENT(IN)   :: coef_in
94    REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay ! pres au milieu de couche (Pa)
95    REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! pression a inter-couche (Pa)
96    REAL, DIMENSION(klon,klev), INTENT(IN)   :: temp  ! temperature
97    REAL, DIMENSION(klon,klev), INTENT(IN)   :: delp
98    REAL, DIMENSION(klon,klev), INTENT(IN)   :: u_old
99    REAL, DIMENSION(klon,klev), INTENT(IN)   :: v_old
100
101! Output arguments
102!****************************************************************************************
103    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef_U_out
104    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef_V_out
105    REAL, DIMENSION(klon), INTENT(OUT)       :: Bcoef_U_out
106    REAL, DIMENSION(klon), INTENT(OUT)       :: Bcoef_V_out
107
108!!! nrlmd le 02/05/2011
109    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Ccoef_U_out
110    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Ccoef_V_out
111    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Dcoef_U_out
112    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Dcoef_V_out
113    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Kcoef_m_out
114    REAL, DIMENSION(klon), INTENT(OUT)       :: alf_1_out
115    REAL, DIMENSION(klon), INTENT(OUT)       :: alf_2_out
116!!!
117
118! Local variables
119!****************************************************************************************
120    REAL, DIMENSION(klon)                    :: u1lay, v1lay
121    INTEGER                                  :: k, i
122!****************************************************************************************
123! Initialize module
124    IF (firstcall) CALL climb_wind_init
125
126!****************************************************************************************
127! Calculate the coefficients C and D in : u(k) = C(k) + D(k)*u(k-1)
128!
129!****************************************************************************************
130! - Define alpha (alf1 and alf2)
131    alf1(:) = 1.0
132    alf2(:) = 1.0 - alf1(:)
133
134! - Calculate the coefficients K
135    Kcoefm(:,:) = 0.0
136    DO k = 2, klev
137       DO i=1,knon
138          Kcoefm(i,k) = coef_in(i,k)*RG*RG*dtime/(pplay(i,k-1)-pplay(i,k)) &
139               *(paprs(i,k)*2/(temp(i,k)+temp(i,k-1))/RD)**2
140       END DO
141    END DO
142
143! - Calculate the coefficients C and D, component "u"
144    CALL calc_coef(knon, Kcoefm(:,:), delp(:,:), &
145         u_old(:,:), alf1(:), alf2(:),  &
146         Ccoef_U(:,:), Dcoef_U(:,:), Acoef_U(:), Bcoef_U(:))
147
148! - Calculate the coefficients C and D, component "v"
149    CALL calc_coef(knon, Kcoefm(:,:), delp(:,:), &
150         v_old(:,:), alf1(:), alf2(:),  &
151         Ccoef_V(:,:), Dcoef_V(:,:), Acoef_V(:), Bcoef_V(:))
152
153!****************************************************************************************
154! 6)
155! Return the first layer in output variables
156!
157!****************************************************************************************
158    Acoef_U_out = Acoef_U
159    Bcoef_U_out = Bcoef_U
160    Acoef_V_out = Acoef_V
161    Bcoef_V_out = Bcoef_V
162
163!****************************************************************************************
164! 7)
165! If Pbl is split, return also the other layers in output variables
166!
167!****************************************************************************************
168!!! jyg le 07/02/2012
169!!jyg       IF (mod(iflag_pbl_split,2) .eq.1) THEN
170       IF (mod(iflag_pbl_split,10) .ge.1) THEN
171!!! nrlmd le 02/05/2011
172    DO k= 1, klev
173      DO i= 1, klon
174        Ccoef_U_out(i,k) = Ccoef_U(i,k)
175        Ccoef_V_out(i,k) = Ccoef_V(i,k)
176        Dcoef_U_out(i,k) = Dcoef_U(i,k)
177        Dcoef_V_out(i,k) = Dcoef_V(i,k)
178        Kcoef_m_out(i,k) = Kcoefm(i,k)
179      ENDDO
180    ENDDO
181    DO i= 1, klon
182      alf_1_out(i)   = alf1(i)
183      alf_2_out(i)   = alf2(i)
184    ENDDO
185!!!     
186       ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
187!!!
188
189  END SUBROUTINE climb_wind_down
190!
191!****************************************************************************************
192!
193  SUBROUTINE calc_coef(knon, Kcoef, delp, X, alfa1, alfa2, Ccoef, Dcoef, Acoef, Bcoef)
194!
195! Find the coefficients C and D in fonction of alfa, K and delp
196USE yomcst_mod_h
197! Input arguments
198!****************************************************************************************
199    INTEGER, INTENT(IN)                      :: knon
200    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Kcoef, delp
201    REAL, DIMENSION(klon,klev), INTENT(IN)   :: X
202    REAL, DIMENSION(klon), INTENT(IN)        :: alfa1, alfa2
203
204! Output arguments
205!****************************************************************************************
206    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef, Bcoef
207    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Ccoef, Dcoef
208 
209! local variables
210!****************************************************************************************
211    INTEGER                                  :: k, i
212    REAL                                     :: buf
213    !****************************************************************************************
214!
215
216! Calculate coefficients C and D at top level, k=klev
217!
218    Ccoef(:,:) = 0.0
219    Dcoef(:,:) = 0.0
220
221    DO i = 1, knon
222       buf = delp(i,klev) + Kcoef(i,klev)
223
224       Ccoef(i,klev) = X(i,klev)*delp(i,klev)/buf
225       Dcoef(i,klev) = Kcoef(i,klev)/buf
226    END DO
227   
228!
229! Calculate coefficients C and D at top level (klev-1) <= k <= 2
230!
231    DO k=(klev-1),2,-1
232       DO i = 1, knon
233          buf = delp(i,k) + Kcoef(i,k) + Kcoef(i,k+1)*(1.-Dcoef(i,k+1))
234         
235          Ccoef(i,k) = (X(i,k)*delp(i,k) + Kcoef(i,k+1)*Ccoef(i,k+1))/buf
236          Dcoef(i,k) = Kcoef(i,k)/buf
237       END DO
238    END DO
239
240!
241! Calculate coeffiecent A and B at surface
242!
243    DO i = 1, knon
244       buf = delp(i,1) + Kcoef(i,2)*(1-Dcoef(i,2))
245       Acoef(i) = (X(i,1)*delp(i,1) + Kcoef(i,2)*Ccoef(i,2))/buf
246       Bcoef(i) = -RG/buf
247    END DO
248
249  END SUBROUTINE calc_coef
250!
251!****************************************************************************************
252!
253
254  SUBROUTINE climb_wind_up(knon, dtime, u_old, v_old, flx_u1, flx_v1,  &
255!!! nrlmd le 02/05/2011
256       Acoef_U_in, Acoef_V_in, Bcoef_U_in, Bcoef_V_in, &
257       Ccoef_U_in, Ccoef_V_in, Dcoef_U_in, Dcoef_V_in, &
258       Kcoef_m_in, &
259!!!
260       flx_u_new, flx_v_new, d_u_new, d_v_new)
261!
262! Diffuse the wind components from the surface layer and up to the top layer.
263! Coefficents A, B, C and D are known from before. Start values for the diffusion are the
264! momentum fluxes at surface.
265!
266! u(k=1) = A + B*flx*dtime
267! u(k)   = C(k) + D(k)*u(k-1)  [2 <= k <= klev]
268!
269!****************************************************************************************
270USE yomcst_mod_h
271USE compbl_mod_h
272! Input arguments
273!****************************************************************************************
274    INTEGER, INTENT(IN)                     :: knon
275    REAL, INTENT(IN)                        :: dtime
276    REAL, DIMENSION(klon,klev), INTENT(IN)  :: u_old
277    REAL, DIMENSION(klon,klev), INTENT(IN)  :: v_old
278    REAL, DIMENSION(klon), INTENT(IN)       :: flx_u1, flx_v1 ! momentum flux
279
280!!! nrlmd le 02/05/2011
281    REAL, DIMENSION(klon), INTENT(IN)       :: Acoef_U_in,Acoef_V_in, Bcoef_U_in, Bcoef_V_in
282    REAL, DIMENSION(klon,klev), INTENT(IN)  :: Ccoef_U_in, Ccoef_V_in, Dcoef_U_in, Dcoef_V_in
283    REAL, DIMENSION(klon,klev), INTENT(IN)  :: Kcoef_m_in
284!!!
285
286! Output arguments
287!****************************************************************************************
288    REAL, DIMENSION(klon,klev), INTENT(OUT) :: flx_u_new, flx_v_new
289    REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_u_new, d_v_new
290
291! Local variables
292!****************************************************************************************
293    REAL, DIMENSION(klon,klev)              :: u_new, v_new
294    INTEGER                                 :: k, i
295!****************************************************************************************
296
297!!! jyg le 07/02/2012
298!!jyg       IF (mod(iflag_pbl_split,2) .eq.1) THEN
299       IF (mod(iflag_pbl_split,10) .ge.1) THEN
300!!! nrlmd le 02/05/2011
301    DO i = 1, knon
302      Acoef_U(i)=Acoef_U_in(i)
303      Acoef_V(i)=Acoef_V_in(i)
304      Bcoef_U(i)=Bcoef_U_in(i)
305      Bcoef_V(i)=Bcoef_V_in(i)
306    ENDDO
307    DO k = 1, klev
308      DO i = 1, knon
309        Ccoef_U(i,k)=Ccoef_U_in(i,k)
310        Ccoef_V(i,k)=Ccoef_V_in(i,k)
311        Dcoef_U(i,k)=Dcoef_U_in(i,k)
312        Dcoef_V(i,k)=Dcoef_V_in(i,k)
313        Kcoefm(i,k)=Kcoef_m_in(i,k)
314      ENDDO
315    ENDDO
316!!!
317       ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
318!!!
319
320! Niveau 1
321    DO i = 1, knon
322       u_new(i,1) = Acoef_U(i) + Bcoef_U(i)*flx_u1(i)*dtime
323       v_new(i,1) = Acoef_V(i) + Bcoef_V(i)*flx_v1(i)*dtime
324    END DO
325
326! Niveau 2 jusqu'au sommet klev
327    DO k = 2, klev
328       DO i=1, knon
329          u_new(i,k) = Ccoef_U(i,k) + Dcoef_U(i,k) * u_new(i,k-1)
330          v_new(i,k) = Ccoef_V(i,k) + Dcoef_V(i,k) * v_new(i,k-1)
331       END DO
332    END DO
333
334!****************************************************************************************
335! Calcul flux
336!
337!== flux_u/v est le flux de moment angulaire (positif vers bas)
338!== dont l'unite est: (kg m/s)/(m**2 s)
339!
340!****************************************************************************************
341!
342    flx_u_new(:,:) = 0.0
343    flx_v_new(:,:) = 0.0
344
345    flx_u_new(1:knon,1)=flx_u1(1:knon)
346    flx_v_new(1:knon,1)=flx_v1(1:knon)
347
348! Niveau 2->klev
349    DO k = 2, klev
350       DO i = 1, knon
351          flx_u_new(i,k) = Kcoefm(i,k)/RG/dtime * &
352               (u_new(i,k)-u_new(i,k-1))
353         
354          flx_v_new(i,k) = Kcoefm(i,k)/RG/dtime * &
355               (v_new(i,k)-v_new(i,k-1))
356       END DO
357    END DO
358
359!****************************************************************************************
360! Calcul tendances
361!
362!****************************************************************************************
363    d_u_new(:,:) = 0.0
364    d_v_new(:,:) = 0.0
365    DO k = 1, klev
366       DO i = 1, knon
367          d_u_new(i,k) = u_new(i,k) - u_old(i,k)
368          d_v_new(i,k) = v_new(i,k) - v_old(i,k)
369       END DO
370    END DO
371
372  END SUBROUTINE climb_wind_up
373!
374!****************************************************************************************
375!
376END MODULE climb_wind_mod
Note: See TracBrowser for help on using the repository browser.