source: LMDZ6/branches/PBLSURF_GPUPORT/libf/phylmd/climb_wind_mod.f90 @ 5914

Last change on this file since 5914 was 5879, checked in by yann meurdesoif, 3 months ago

GPU port of climb_wind_up + climb_wind_down

YM
  • 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.3 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, climb_wind_init
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!****************************************************************************************
42! Initialize module
43   
44    IF (firstcall) THEN
45
46      !****************************************************************************************
47      ! Allocation of global module variables
48      !
49      !****************************************************************************************
50 
51      ALLOCATE(alf1(klon), stat=ierr)
52      IF (ierr /= 0) CALL abort_physic(modname,'Pb in allocate alf1',1)
53      alf1(:) = 0.
54 
55      ALLOCATE(alf2(klon), stat=ierr)
56      IF (ierr /= 0) CALL abort_physic(modname,'Pb in allocate alf2',1)
57      alf2(:) = 0.
58 
59      ALLOCATE(Kcoefm(klon,klev), stat=ierr)
60      IF (ierr /= 0) CALL abort_physic(modname,'Pb in allocate Kcoefm',1)
61      Kcoefm(:,:) = 0.
62 
63      ALLOCATE(Ccoef_U(klon,klev), stat=ierr)
64      IF (ierr /= 0) CALL abort_physic(modname,'Pb in allocate Ccoef_U',1)
65      Ccoef_U(:,:) = 0.
66 
67      ALLOCATE(Dcoef_U(klon,klev), stat=ierr)
68      IF (ierr /= 0) CALL abort_physic(modname,'Pb in allocation Dcoef_U',1)
69      Dcoef_U(:,:) = 0.
70 
71      ALLOCATE(Ccoef_V(klon,klev), stat=ierr)
72      IF (ierr /= 0) CALL abort_physic(modname,'Pb in allocation Ccoef_V',1)
73      Ccoef_V(:,:) = 0.
74 
75      ALLOCATE(Dcoef_V(klon,klev), stat=ierr)
76      IF (ierr /= 0) CALL abort_physic(modname,'Pb in allocation Dcoef_V',1)
77      Dcoef_V(:,:) = 0.
78 
79      ALLOCATE(Acoef_U(klon), Bcoef_U(klon), Acoef_V(klon), Bcoef_V(klon), STAT=ierr)
80      IF ( ierr /= 0 )  PRINT*,' pb in allloc Acoef_U and Bcoef_U, ierr=', ierr
81      Acoef_U(:) = 0. ; Bcoef_U(:) = 0. ; Acoef_V(:) = 0. ; Bcoef_V(:) = 0. ;
82 
83      firstcall=.FALSE.
84   
85    ENDIF
86
87  END SUBROUTINE climb_wind_init
88!
89!****************************************************************************************
90!
91  SUBROUTINE climb_wind_down(knon, ni, dtime, coef_in, pplay, paprs, temp, delp, u_old, v_old, &
92!!! nrlmd le 02/05/2011
93       Ccoef_U_out, Ccoef_V_out, Dcoef_U_out, Dcoef_V_out, &
94       Kcoef_m_out, alf_1_out, alf_2_out, &
95!!!
96       Acoef_U_out, Acoef_V_out, Bcoef_U_out, Bcoef_V_out)
97!$gpum horizontal knon
98
99!
100! This routine calculates for the wind components u and v,
101! recursivly the coefficients C and D in equation
102! X(k) = C(k) + D(k)*X(k-1), X=[u,v], k=[1,klev] is the vertical layer.
103!
104!
105USE yomcst_mod_h
106USE compbl_mod_h
107! Input arguments
108!****************************************************************************************
109    INTEGER, INTENT(IN)                      :: knon
110    INTEGER, INTENT(IN)                      :: ni(knon)
111    REAL, INTENT(IN)                         :: dtime
112    REAL, DIMENSION(knon,klev), INTENT(IN)   :: coef_in
113    REAL, DIMENSION(knon,klev), INTENT(IN)   :: pplay ! pres au milieu de couche (Pa)
114    REAL, DIMENSION(knon,klev+1), INTENT(IN) :: paprs ! pression a inter-couche (Pa)
115    REAL, DIMENSION(knon,klev), INTENT(IN)   :: temp  ! temperature
116    REAL, DIMENSION(knon,klev), INTENT(IN)   :: delp
117    REAL, DIMENSION(knon,klev), INTENT(IN)   :: u_old
118    REAL, DIMENSION(knon,klev), INTENT(IN)   :: v_old
119
120! Output arguments
121!****************************************************************************************
122    REAL, DIMENSION(knon), INTENT(OUT)       :: Acoef_U_out
123    REAL, DIMENSION(knon), INTENT(OUT)       :: Acoef_V_out
124    REAL, DIMENSION(knon), INTENT(OUT)       :: Bcoef_U_out
125    REAL, DIMENSION(knon), INTENT(OUT)       :: Bcoef_V_out
126
127!!! nrlmd le 02/05/2011
128    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: Ccoef_U_out
129    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: Ccoef_V_out
130    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: Dcoef_U_out
131    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: Dcoef_V_out
132    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: Kcoef_m_out
133    REAL, DIMENSION(knon), INTENT(OUT)       :: alf_1_out
134    REAL, DIMENSION(knon), INTENT(OUT)       :: alf_2_out
135!!!
136
137! Local variables
138!****************************************************************************************
139    REAL :: yalf1(knon)
140    REAL :: yalf2(knon)
141    REAL :: yKcoefm(knon,klev)
142    REAL :: yCcoef_U(knon,klev)
143    REAL :: yDcoef_U(knon,klev)
144    REAL :: yCcoef_V(knon,klev)
145    REAL :: yDcoef_V(knon,klev)
146    REAL :: yAcoef_U(knon), yBcoef_U(knon), yAcoef_V(knon), yBcoef_V(knon)
147
148    REAL, DIMENSION(knon)                    :: u1lay, v1lay
149    INTEGER                                  :: k, i, j
150
151
152!****************************************************************************************
153! Calculate the coefficients C and D in : u(k) = C(k) + D(k)*u(k-1)
154!
155!****************************************************************************************
156! - Define alpha (alf1 and alf2)
157    yalf1(:) = 1.0
158    yalf2(:) = 1.0 - yalf1(:)
159
160! - Calculate the coefficients K
161    yKcoefm(:,:) = 0.0
162    DO k = 2, klev
163       DO i=1,knon
164          yKcoefm(i,k) = coef_in(i,k)*RG*RG*dtime/(pplay(i,k-1)-pplay(i,k)) &
165               *(paprs(i,k)*2/(temp(i,k)+temp(i,k-1))/RD)**2
166       END DO
167    END DO
168
169! - Calculate the coefficients C and D, component "u"
170    CALL calc_coef(knon, yKcoefm, delp, &
171         u_old, yalf1, yalf2,  &
172         yCcoef_U, yDcoef_U, yAcoef_U, yBcoef_U)
173
174! - Calculate the coefficients C and D, component "v"
175    CALL calc_coef(knon, yKcoefm, delp, &
176         v_old, yalf1, yalf2,  &
177         yCcoef_V, yDcoef_V, yAcoef_V, yBcoef_V)
178
179!****************************************************************************************
180! 6)
181! Return the first layer in output variables
182!
183!****************************************************************************************
184    Acoef_U_out = yAcoef_U
185    Bcoef_U_out = yBcoef_U
186    Acoef_V_out = yAcoef_V
187    Bcoef_V_out = yBcoef_V
188
189!****************************************************************************************
190! 7)
191! If Pbl is split, return also the other layers in output variables
192!
193!****************************************************************************************
194!!! jyg le 07/02/2012
195!!jyg       IF (mod(iflag_pbl_split,2) .eq.1) THEN
196       IF (mod(iflag_pbl_split,10) .ge.1) THEN
197!!! nrlmd le 02/05/2011
198    DO k= 1, klev
199      DO i= 1, knon
200        Ccoef_U_out(i,k) = yCcoef_U(i,k)
201        Ccoef_V_out(i,k) = yCcoef_V(i,k)
202        Dcoef_U_out(i,k) = yDcoef_U(i,k)
203        Dcoef_V_out(i,k) = yDcoef_V(i,k)
204        Kcoef_m_out(i,k) = yKcoefm(i,k)
205      ENDDO
206    ENDDO
207    DO i= 1, knon
208      alf_1_out(i)   = yalf1(i)
209      alf_2_out(i)   = yalf2(i)
210    ENDDO
211!!!     
212       ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
213!!!
214    DO j= 1, knon
215      i=ni(j)
216      Acoef_U(i) = yAcoef_U(j)
217      Bcoef_U(i) = yBcoef_U(j)
218      Acoef_V(i) = yAcoef_V(j)
219      Bcoef_V(i) = yBcoef_V(j)
220    ENDDO
221
222    DO k= 1, klev
223      DO j= 1, knon
224        i=ni(j)
225        Ccoef_U(i,k) = yCcoef_U(j,k)
226        Ccoef_V(i,k) = yCcoef_V(j,k)
227        Dcoef_U(i,k) = yDcoef_U(j,k)
228        Dcoef_V(i,k) = yDcoef_V(j,k)
229        Kcoefm(i,k) = yKcoefm(j,k)
230        alf1(i)   = yalf1(j)
231        alf2(i)   = yalf2(j)
232      ENDDO
233    ENDDO
234  END SUBROUTINE climb_wind_down
235!
236!****************************************************************************************
237!
238  SUBROUTINE calc_coef(knon, Kcoef, delp, X, alfa1, alfa2, Ccoef, Dcoef, Acoef, Bcoef)
239!$gpum horizontal knon
240
241!
242! Find the coefficients C and D in fonction of alfa, K and delp
243USE yomcst_mod_h
244! Input arguments
245!****************************************************************************************
246    INTEGER, INTENT(IN)                      :: knon
247    REAL, DIMENSION(knon,klev), INTENT(IN)   :: Kcoef, delp
248    REAL, DIMENSION(knon,klev), INTENT(IN)   :: X
249    REAL, DIMENSION(knon), INTENT(IN)        :: alfa1, alfa2
250
251! Output arguments
252!****************************************************************************************
253    REAL, DIMENSION(knon), INTENT(OUT)       :: Acoef, Bcoef
254    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: Ccoef, Dcoef
255 
256! local variables
257!****************************************************************************************
258    INTEGER                                  :: k, i, j
259    REAL                                     :: buf
260    !****************************************************************************************
261!
262
263! Calculate coefficients C and D at top level, k=klev
264!
265    Ccoef(:,:) = 0.0
266    Dcoef(:,:) = 0.0
267
268    DO i = 1, knon
269       buf = delp(i,klev) + Kcoef(i,klev)
270
271       Ccoef(i,klev) = X(i,klev)*delp(i,klev)/buf
272       Dcoef(i,klev) = Kcoef(i,klev)/buf
273    END DO
274   
275!
276! Calculate coefficients C and D at top level (klev-1) <= k <= 2
277!
278    DO k=(klev-1),2,-1
279       DO i = 1, knon
280          buf = delp(i,k) + Kcoef(i,k) + Kcoef(i,k+1)*(1.-Dcoef(i,k+1))
281         
282          Ccoef(i,k) = (X(i,k)*delp(i,k) + Kcoef(i,k+1)*Ccoef(i,k+1))/buf
283          Dcoef(i,k) = Kcoef(i,k)/buf
284       END DO
285    END DO
286
287!
288! Calculate coeffiecent A and B at surface
289!
290    DO i = 1, knon
291       buf = delp(i,1) + Kcoef(i,2)*(1-Dcoef(i,2))
292       Acoef(i) = (X(i,1)*delp(i,1) + Kcoef(i,2)*Ccoef(i,2))/buf
293       Bcoef(i) = -RG/buf
294    END DO
295
296  END SUBROUTINE calc_coef
297!
298!****************************************************************************************
299!
300
301  SUBROUTINE climb_wind_up(knon, ni, dtime, u_old, v_old, flx_u1, flx_v1,  &
302!!! nrlmd le 02/05/2011
303       Acoef_U_in, Acoef_V_in, Bcoef_U_in, Bcoef_V_in, &
304       Ccoef_U_in, Ccoef_V_in, Dcoef_U_in, Dcoef_V_in, &
305       Kcoef_m_in, &
306!!!
307       flx_u_new, flx_v_new, d_u_new, d_v_new)
308!$gpum horizontal knon
309
310!
311! Diffuse the wind components from the surface layer and up to the top layer.
312! Coefficents A, B, C and D are known from before. Start values for the diffusion are the
313! momentum fluxes at surface.
314!
315! u(k=1) = A + B*flx*dtime
316! u(k)   = C(k) + D(k)*u(k-1)  [2 <= k <= klev]
317!
318!****************************************************************************************
319USE yomcst_mod_h
320USE compbl_mod_h
321! Input arguments
322!****************************************************************************************
323    INTEGER, INTENT(IN)                     :: knon
324    INTEGER, INTENT(IN)                     :: ni(knon)
325    REAL, INTENT(IN)                        :: dtime
326    REAL, DIMENSION(knon,klev), INTENT(IN)  :: u_old
327    REAL, DIMENSION(knon,klev), INTENT(IN)  :: v_old
328    REAL, DIMENSION(knon), INTENT(IN)       :: flx_u1, flx_v1 ! momentum flux
329
330!!! nrlmd le 02/05/2011
331    REAL, DIMENSION(knon), INTENT(IN)       :: Acoef_U_in,Acoef_V_in, Bcoef_U_in, Bcoef_V_in
332    REAL, DIMENSION(knon,klev), INTENT(IN)  :: Ccoef_U_in, Ccoef_V_in, Dcoef_U_in, Dcoef_V_in
333    REAL, DIMENSION(knon,klev), INTENT(IN)  :: Kcoef_m_in
334!!!
335
336! Output arguments
337!****************************************************************************************
338    REAL, DIMENSION(knon,klev), INTENT(OUT) :: flx_u_new, flx_v_new
339    REAL, DIMENSION(knon,klev), INTENT(OUT) :: d_u_new, d_v_new
340
341! Local variables
342!****************************************************************************************
343    REAL :: yalf1(knon)
344    REAL :: yalf2(knon)
345    REAL :: yKcoefm(knon,klev)
346    REAL :: yCcoef_U(knon,klev)
347    REAL :: yDcoef_U(knon,klev)
348    REAL :: yCcoef_V(knon,klev)
349    REAL :: yDcoef_V(knon,klev)
350    REAL :: yAcoef_U(knon), yBcoef_U(knon), yAcoef_V(knon), yBcoef_V(knon)
351
352    REAL, DIMENSION(knon,klev)              :: u_new, v_new
353    INTEGER                                 :: k, i, j
354!****************************************************************************************
355
356!!! jyg le 07/02/2012
357!!jyg       IF (mod(iflag_pbl_split,2) .eq.1) THEN
358       IF (mod(iflag_pbl_split,10) .ge.1) THEN
359!!! nrlmd le 02/05/2011
360    DO i = 1, knon
361      yAcoef_U(i)=Acoef_U_in(i)
362      yAcoef_V(i)=Acoef_V_in(i)
363      yBcoef_U(i)=Bcoef_U_in(i)
364      yBcoef_V(i)=Bcoef_V_in(i)
365    ENDDO
366    DO k = 1, klev
367      DO i = 1, knon
368        yCcoef_U(i,k)=Ccoef_U_in(i,k)
369        yCcoef_V(i,k)=Ccoef_V_in(i,k)
370        yDcoef_U(i,k)=Dcoef_U_in(i,k)
371        yDcoef_V(i,k)=Dcoef_V_in(i,k)
372        yKcoefm(i,k)=Kcoef_m_in(i,k)
373      ENDDO
374    ENDDO
375   ELSE
376    DO j = 1, knon
377      i=ni(j)
378      yAcoef_U(j)=Acoef_U(i)
379      yAcoef_V(j)=Acoef_V(i)
380      yBcoef_U(j)=Bcoef_U(i)
381      yBcoef_V(j)=Bcoef_V(i)
382    ENDDO
383    DO k = 1, klev
384      DO j = 1, knon
385        i=ni(j)
386        yCcoef_U(j,k)=Ccoef_U(i,k)
387        yCcoef_V(j,k)=Ccoef_V(i,k)
388        yDcoef_U(j,k)=Dcoef_U(i,k)
389        yDcoef_V(j,k)=Dcoef_V(i,k)
390        yKcoefm(j,k)=Kcoefm(i,k)
391      ENDDO
392    ENDDO   
393!!!
394       ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
395!!!
396
397! Niveau 1
398    DO i = 1, knon
399       u_new(i,1) = yAcoef_U(i) + yBcoef_U(i)*flx_u1(i)*dtime
400       v_new(i,1) = yAcoef_V(i) + yBcoef_V(i)*flx_v1(i)*dtime
401    END DO
402
403! Niveau 2 jusqu'au sommet klev
404    DO k = 2, klev
405       DO i=1, knon
406          u_new(i,k) = yCcoef_U(i,k) + yDcoef_U(i,k) * u_new(i,k-1)
407          v_new(i,k) = yCcoef_V(i,k) + yDcoef_V(i,k) * v_new(i,k-1)
408       END DO
409    END DO
410
411!****************************************************************************************
412! Calcul flux
413!
414!== flux_u/v est le flux de moment angulaire (positif vers bas)
415!== dont l'unite est: (kg m/s)/(m**2 s)
416!
417!****************************************************************************************
418!
419    flx_u_new(:,:) = 0.0
420    flx_v_new(:,:) = 0.0
421
422    flx_u_new(1:knon,1)=flx_u1(1:knon)
423    flx_v_new(1:knon,1)=flx_v1(1:knon)
424
425! Niveau 2->klev
426    DO k = 2, klev
427       DO i = 1, knon
428          flx_u_new(i,k) = yKcoefm(i,k)/RG/dtime * &
429               (u_new(i,k)-u_new(i,k-1))
430         
431          flx_v_new(i,k) = yKcoefm(i,k)/RG/dtime * &
432               (v_new(i,k)-v_new(i,k-1))
433       END DO
434    END DO
435
436!****************************************************************************************
437! Calcul tendances
438!
439!****************************************************************************************
440    d_u_new(:,:) = 0.0
441    d_v_new(:,:) = 0.0
442    DO k = 1, klev
443       DO i = 1, knon
444          d_u_new(i,k) = u_new(i,k) - u_old(i,k)
445          d_v_new(i,k) = v_new(i,k) - v_old(i,k)
446       END DO
447    END DO
448
449    DO j = 1, knon
450      i=ni(j)
451      Acoef_U(i)=yAcoef_U(j)
452      Acoef_V(i)=yAcoef_V(j)
453      Bcoef_U(i)=yBcoef_U(j)
454      Bcoef_V(i)=yBcoef_V(j)
455    ENDDO
456
457    DO k = 1, klev
458      DO j = 1, knon
459        i=ni(j)
460        Ccoef_U(i,k) = yCcoef_U(j,k)
461        Ccoef_V(i,k) = yCcoef_V(j,k)
462        Dcoef_U(i,k) = yDcoef_U(j,k)
463        Dcoef_V(i,k) = yDcoef_V(j,k)
464        Kcoefm(i,k) = yKcoefm(j,k)
465      ENDDO
466    ENDDO   
467 
468  END SUBROUTINE climb_wind_up
469!
470!****************************************************************************************
471!
472END MODULE climb_wind_mod
Note: See TracBrowser for help on using the repository browser.