source: codes/icosagcm/trunk/src/caldyn_kernels_base.f90 @ 366

Last change on this file since 366 was 366, checked in by dubos, 9 years ago

Progress towards NH
File size: 8.9 KB
Line 
1MODULE caldyn_kernels_base_mod
2  USE icosa
3  USE transfert_mod
4  IMPLICIT NONE
5  PRIVATE
6
7  LOGICAL, PARAMETER, PUBLIC :: DEC = .TRUE.
8
9  INTEGER, PARAMETER,PUBLIC :: energy=1, enstrophy=2
10  TYPE(t_field),POINTER,PUBLIC :: f_out_u(:), f_qu(:), f_qv(:)
11  REAL(rstd),SAVE,POINTER :: out_u(:,:), p(:,:), qu(:,:)
12  !$OMP THREADPRIVATE(out_u, p, qu)
13
14  ! temporary shared variables for caldyn
15  TYPE(t_field),POINTER,PUBLIC :: f_pk(:),f_wwuu(:),f_planetvel(:)
16
17  INTEGER, PUBLIC :: caldyn_conserv
18  !$OMP THREADPRIVATE(caldyn_conserv)
19
20  TYPE(t_message),PUBLIC :: req_ps, req_mass, req_theta_rhodz, req_u, req_qu, req_geopot, req_w
21
22  PUBLIC :: compute_geopot, compute_caldyn_vert
23
24CONTAINS
25
26  !**************************** Geopotential *****************************
27
28  SUBROUTINE compute_geopot(ps,rhodz,theta, pk,geopot) 
29    USE icosa
30    USE disvert_mod
31    USE exner_mod
32    USE trace
33    USE omp_para
34    IMPLICIT NONE
35    REAL(rstd),INTENT(INOUT) :: ps(iim*jjm)
36    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
37    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm)    ! potential temperature
38    REAL(rstd),INTENT(OUT)   :: pk(iim*jjm,llm)       ! Exner function (compressible) /Lagrange multiplier (Boussinesq)
39    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) ! geopotential
40
41    INTEGER :: i,j,ij,l
42    REAL(rstd) :: p_ik, exner_ik
43    INTEGER    :: ij_omp_begin_ext, ij_omp_end_ext
44
45
46    CALL trace_start("compute_geopot")
47
48    CALL distrib_level(ij_end_ext-ij_begin_ext+1,ij_omp_begin_ext,ij_omp_end_ext)
49    ij_omp_begin_ext=ij_omp_begin_ext+ij_begin_ext-1
50    ij_omp_end_ext=ij_omp_end_ext+ij_begin_ext-1
51
52    IF(caldyn_eta==eta_mass .AND. .NOT. DEC) THEN
53!!! Compute exner function and geopotential
54       DO l = 1,llm
55          !DIR$ SIMD
56          DO ij=ij_omp_begin_ext,ij_omp_end_ext
57             p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later
58             exner_ik = cpp * (p_ik/preff) ** kappa
59             pk(ij,l) = exner_ik
60             ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
61             geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik
62          ENDDO
63       ENDDO
64    ELSE 
65       ! We are using DEC or a Lagrangian vertical coordinate
66       ! Pressure is computed first top-down (temporarily stored in pk)
67       ! Then Exner pressure and geopotential are computed bottom-up
68       ! Works also when caldyn_eta=eta_mass         
69
70       IF(boussinesq) THEN ! compute geopotential and pk=Lagrange multiplier
71          ! specific volume 1 = dphi/g/rhodz
72          !         IF (is_omp_level_master) THEN ! no openMP on vertical due to dependency
73          DO l = 1,llm
74             !DIR$ SIMD
75             DO ij=ij_omp_begin_ext,ij_omp_end_ext         
76                geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l)
77             ENDDO
78          ENDDO
79          ! use hydrostatic balance with theta*rhodz to find pk (Lagrange multiplier=pressure)
80          ! uppermost layer
81          !DIR$ SIMD
82          DO ij=ij_begin_ext,ij_end_ext         
83             pk(ij,llm) = ptop + (.5*g)*theta(ij,llm)*rhodz(ij,llm)
84          END DO
85          ! other layers
86          DO l = llm-1, 1, -1
87             !          !$OMP DO SCHEDULE(STATIC)
88             !DIR$ SIMD
89             DO ij=ij_begin_ext,ij_end_ext         
90                pk(ij,l) = pk(ij,l+1) + (.5*g)*(theta(ij,l)*rhodz(ij,l)+theta(ij,l+1)*rhodz(ij,l+1))
91             END DO
92          END DO
93          ! now pk contains the Lagrange multiplier (pressure)
94       ELSE ! non-Boussinesq, compute geopotential and Exner pressure
95          ! uppermost layer
96
97          !DIR$ SIMD
98          DO ij=ij_omp_begin_ext,ij_omp_end_ext
99             pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm)
100          END DO
101          ! other layers
102          DO l = llm-1, 1, -1
103             !DIR$ SIMD
104             DO ij=ij_omp_begin_ext,ij_omp_end_ext         
105                pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1))
106             END DO
107          END DO
108          ! surface pressure (for diagnostics)
109          IF(caldyn_eta==eta_lag) THEN
110             DO ij=ij_omp_begin_ext,ij_omp_end_ext         
111                ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1)
112             END DO
113          END IF
114          ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
115          DO l = 1,llm
116             !DIR$ SIMD
117             DO ij=ij_omp_begin_ext,ij_omp_end_ext         
118                p_ik = pk(ij,l)
119                exner_ik = cpp * (p_ik/preff) ** kappa
120                pk(ij,l) = exner_ik
121                geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 
122             ENDDO
123          ENDDO
124       END IF
125
126    END IF
127
128    !ym flush geopot
129    !$OMP BARRIER
130
131    CALL trace_end("compute_geopot")
132
133  END SUBROUTINE compute_geopot
134
135  SUBROUTINE compute_caldyn_vert(u,theta,rhodz,convm, wflux,wwuu, dps,dtheta_rhodz,du)
136    USE icosa
137    USE disvert_mod
138    USE exner_mod
139    USE trace
140    USE omp_para
141    IMPLICIT NONE
142    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
143    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)
144    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
145    REAL(rstd),INTENT(INOUT)  :: convm(iim*jjm,llm)  ! mass flux convergence
146
147    REAL(rstd),INTENT(INOUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s)
148    REAL(rstd),INTENT(INOUT) :: wwuu(iim*3*jjm,llm+1)
149    REAL(rstd),INTENT(INOUT) :: du(iim*3*jjm,llm)
150    REAL(rstd),INTENT(INOUT) :: dtheta_rhodz(iim*jjm,llm)
151    REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
152
153    ! temporary variable   
154    INTEGER :: i,j,ij,l
155    REAL(rstd) :: p_ik, exner_ik
156    INTEGER    :: ij_omp_begin, ij_omp_end
157
158    CALL trace_start("compute_caldyn_vert")
159
160    CALL distrib_level(ij_end-ij_begin+1,ij_omp_begin,ij_omp_end)
161    ij_omp_begin=ij_omp_begin+ij_begin-1
162    ij_omp_end=ij_omp_end+ij_begin-1
163
164    !    REAL(rstd) :: wwuu(iim*3*jjm,llm+1) ! tmp var, don't know why but gain 30% on the whole code in opemp
165    ! need to be understood
166
167    !    wwuu=wwuu_out
168    CALL trace_start("compute_caldyn_vert")
169
170    !$OMP BARRIER   
171!!! cumulate mass flux convergence from top to bottom
172    !  IF (is_omp_level_master) THEN
173    DO  l = llm-1, 1, -1
174       !    IF (caldyn_conserv==energy) CALL test_message(req_qu)
175
176!!$OMP DO SCHEDULE(STATIC)
177       !DIR$ SIMD
178       DO ij=ij_omp_begin,ij_omp_end         
179          convm(ij,l) = convm(ij,l) + convm(ij,l+1)
180       ENDDO
181    ENDDO
182    !  ENDIF
183
184    !$OMP BARRIER
185    ! FLUSH on convm
186!!!!!!!!!!!!!!!!!!!!!!!!!
187
188    ! compute dps
189    IF (is_omp_first_level) THEN
190       !DIR$ SIMD
191       DO ij=ij_begin,ij_end         
192          ! dps/dt = -int(div flux)dz
193          IF(DEC) THEN
194             dps(ij) = convm(ij,1)
195          ELSE
196             dps(ij) = convm(ij,1) * g 
197          END IF
198       ENDDO
199    ENDIF
200
201!!! Compute vertical mass flux (l=1,llm+1 done by caldyn_BC)
202    DO l=ll_beginp1,ll_end
203       !    IF (caldyn_conserv==energy) CALL test_message(req_qu)
204       !DIR$ SIMD
205       DO ij=ij_begin,ij_end         
206          ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
207          ! => w>0 for upward transport
208          wflux( ij, l ) = bp(l) * convm( ij, 1 ) - convm( ij, l ) 
209       ENDDO
210    ENDDO
211
212    !--> flush wflux 
213    !$OMP BARRIER
214
215    DO l=ll_begin,ll_endm1
216       !DIR$ SIMD
217       DO ij=ij_begin,ij_end         
218          dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -   0.5 * (  wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1)))
219       ENDDO
220    ENDDO
221
222    DO l=ll_beginp1,ll_end
223       !DIR$ SIMD
224       DO ij=ij_begin,ij_end         
225          dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  +   0.5 * ( wflux(ij,l  ) * (theta(ij,l-1) + theta(ij,l) ) )
226       ENDDO
227    ENDDO
228
229
230    ! Compute vertical transport
231    DO l=ll_beginp1,ll_end
232       !DIR$ SIMD
233       DO ij=ij_begin,ij_end         
234          wwuu(ij+u_right,l) = 0.5*( wflux(ij,l) + wflux(ij+t_right,l)) * (u(ij+u_right,l) - u(ij+u_right,l-1))
235          wwuu(ij+u_lup,l) = 0.5* ( wflux(ij,l) + wflux(ij+t_lup,l)) * (u(ij+u_lup,l) - u(ij+u_lup,l-1))
236          wwuu(ij+u_ldown,l) = 0.5*( wflux(ij,l) + wflux(ij+t_ldown,l)) * (u(ij+u_ldown,l) - u(ij+u_ldown,l-1))
237       ENDDO
238    ENDDO
239
240    !--> flush wwuu 
241    !$OMP BARRIER
242
243    ! Add vertical transport to du
244    DO l=ll_begin,ll_end
245       !DIR$ SIMD
246       DO ij=ij_begin,ij_end         
247          du(ij+u_right, l )   = du(ij+u_right,l)  - (wwuu(ij+u_right,l+1)+ wwuu(ij+u_right,l)) / (rhodz(ij,l)+rhodz(ij+t_right,l))
248          du(ij+u_lup, l )     = du(ij+u_lup,l)    - (wwuu(ij+u_lup,l+1)  + wwuu(ij+u_lup,l))   / (rhodz(ij,l)+rhodz(ij+t_lup,l))
249          du(ij+u_ldown, l )   = du(ij+u_ldown,l)  - (wwuu(ij+u_ldown,l+1)+ wwuu(ij+u_ldown,l)) / (rhodz(ij,l)+rhodz(ij+t_ldown,l))
250       ENDDO
251    ENDDO
252
253    !  DO l=ll_beginp1,ll_end
254    !!DIR$ SIMD
255    !      DO ij=ij_begin,ij_end         
256    !        wwuu_out(ij+u_right,l) = wwuu(ij+u_right,l)
257    !        wwuu_out(ij+u_lup,l)   = wwuu(ij+u_lup,l)
258    !        wwuu_out(ij+u_ldown,l) = wwuu(ij+u_ldown,l)
259    !     ENDDO
260    !   ENDDO
261
262    CALL trace_end("compute_caldyn_vert")
263
264  END SUBROUTINE compute_caldyn_vert
265
266END MODULE caldyn_kernels_base_mod
Note: See TracBrowser for help on using the repository browser.