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

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

Introduced DEC convention for velocity - HEVI only - cleanup to follow
File size: 9.0 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
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
54!!! Compute exner function and geopotential
55       DO l = 1,llm
56          !DIR$ SIMD
57          DO ij=ij_omp_begin_ext,ij_omp_end_ext
58             p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later
59             exner_ik = cpp * (p_ik/preff) ** kappa
60             pk(ij,l) = exner_ik
61             ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
62             geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik
63          ENDDO
64       ENDDO
65    ELSE 
66       ! We are using a Lagrangian vertical coordinate
67       ! Pressure must be computed first top-down (temporarily stored in pk)
68       ! Then Exner pressure and geopotential are computed bottom-up
69       ! Notice that the computation below should work also when caldyn_eta=eta_mass         
70
71       IF(boussinesq) THEN ! compute geopotential and pk=Lagrange multiplier
72          ! specific volume 1 = dphi/g/rhodz
73          !         IF (is_omp_level_master) THEN ! no openMP on vertical due to dependency
74          DO l = 1,llm
75             !DIR$ SIMD
76             DO ij=ij_omp_begin_ext,ij_omp_end_ext         
77                geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l)
78             ENDDO
79          ENDDO
80          ! use hydrostatic balance with theta*rhodz to find pk (Lagrange multiplier=pressure)
81          ! uppermost layer
82          !DIR$ SIMD
83          DO ij=ij_begin_ext,ij_end_ext         
84             pk(ij,llm) = ptop + (.5*g)*theta(ij,llm)*rhodz(ij,llm)
85          END DO
86          ! other layers
87          DO l = llm-1, 1, -1
88             !          !$OMP DO SCHEDULE(STATIC)
89             !DIR$ SIMD
90             DO ij=ij_begin_ext,ij_end_ext         
91                pk(ij,l) = pk(ij,l+1) + (.5*g)*(theta(ij,l)*rhodz(ij,l)+theta(ij,l+1)*rhodz(ij,l+1))
92             END DO
93          END DO
94          ! now pk contains the Lagrange multiplier (pressure)
95       ELSE ! non-Boussinesq, compute geopotential and Exner pressure
96          ! uppermost layer
97
98          !DIR$ SIMD
99          DO ij=ij_omp_begin_ext,ij_omp_end_ext
100             pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm)
101          END DO
102          ! other layers
103          DO l = llm-1, 1, -1
104             !DIR$ SIMD
105             DO ij=ij_omp_begin_ext,ij_omp_end_ext         
106                pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1))
107             END DO
108          END DO
109          ! surface pressure (for diagnostics)
110          IF(caldyn_eta==eta_lag) THEN
111             DO ij=ij_omp_begin_ext,ij_omp_end_ext         
112                ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1)
113             END DO
114          END IF
115          ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
116          DO l = 1,llm
117             !DIR$ SIMD
118             DO ij=ij_omp_begin_ext,ij_omp_end_ext         
119                p_ik = pk(ij,l)
120                exner_ik = cpp * (p_ik/preff) ** kappa
121                pk(ij,l) = exner_ik
122                geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 
123             ENDDO
124          ENDDO
125       END IF
126
127    END IF
128
129    !ym flush geopot
130    !$OMP BARRIER
131
132    CALL trace_end("compute_geopot")
133
134  END SUBROUTINE compute_geopot
135
136  SUBROUTINE compute_caldyn_vert(u,theta,rhodz,convm, wflux,wwuu, dps,dtheta_rhodz,du)
137    USE icosa
138    USE disvert_mod
139    USE exner_mod
140    USE trace
141    USE omp_para
142    IMPLICIT NONE
143    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
144    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)
145    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
146    REAL(rstd),INTENT(INOUT)  :: convm(iim*jjm,llm)  ! mass flux convergence
147
148    REAL(rstd),INTENT(INOUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s)
149    REAL(rstd),INTENT(INOUT) :: wwuu(iim*3*jjm,llm+1)
150    REAL(rstd),INTENT(INOUT) :: du(iim*3*jjm,llm)
151    REAL(rstd),INTENT(INOUT) :: dtheta_rhodz(iim*jjm,llm)
152    REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
153
154    ! temporary variable   
155    INTEGER :: i,j,ij,l
156    REAL(rstd) :: p_ik, exner_ik
157    INTEGER    :: ij_omp_begin, ij_omp_end
158
159    CALL trace_start("compute_caldyn_vert")
160
161    CALL distrib_level(ij_end-ij_begin+1,ij_omp_begin,ij_omp_end)
162    ij_omp_begin=ij_omp_begin+ij_begin-1
163    ij_omp_end=ij_omp_end+ij_begin-1
164
165    !    REAL(rstd) :: wwuu(iim*3*jjm,llm+1) ! tmp var, don't know why but gain 30% on the whole code in opemp
166    ! need to be understood
167
168    !    wwuu=wwuu_out
169    CALL trace_start("compute_caldyn_vert")
170
171    !$OMP BARRIER   
172!!! cumulate mass flux convergence from top to bottom
173    !  IF (is_omp_level_master) THEN
174    DO  l = llm-1, 1, -1
175       !    IF (caldyn_conserv==energy) CALL test_message(req_qu)
176
177!!$OMP DO SCHEDULE(STATIC)
178       !DIR$ SIMD
179       DO ij=ij_omp_begin,ij_omp_end         
180          convm(ij,l) = convm(ij,l) + convm(ij,l+1)
181       ENDDO
182    ENDDO
183    !  ENDIF
184
185    !$OMP BARRIER
186    ! FLUSH on convm
187!!!!!!!!!!!!!!!!!!!!!!!!!
188
189    ! compute dps
190    IF (is_omp_first_level) THEN
191       !DIR$ SIMD
192       DO ij=ij_begin,ij_end         
193          ! dps/dt = -int(div flux)dz
194          IF(DEC) THEN
195             dps(ij) = convm(ij,1)
196          ELSE
197             dps(ij) = convm(ij,1) * g 
198          END IF
199       ENDDO
200    ENDIF
201
202!!! Compute vertical mass flux (l=1,llm+1 done by caldyn_BC)
203    DO l=ll_beginp1,ll_end
204       !    IF (caldyn_conserv==energy) CALL test_message(req_qu)
205       !DIR$ SIMD
206       DO ij=ij_begin,ij_end         
207          ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
208          ! => w>0 for upward transport
209          wflux( ij, l ) = bp(l) * convm( ij, 1 ) - convm( ij, l ) 
210       ENDDO
211    ENDDO
212
213    !--> flush wflux 
214    !$OMP BARRIER
215
216    DO l=ll_begin,ll_endm1
217       !DIR$ SIMD
218       DO ij=ij_begin,ij_end         
219          dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -   0.5 * (  wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1)))
220       ENDDO
221    ENDDO
222
223    DO l=ll_beginp1,ll_end
224       !DIR$ SIMD
225       DO ij=ij_begin,ij_end         
226          dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  +   0.5 * ( wflux(ij,l  ) * (theta(ij,l-1) + theta(ij,l) ) )
227       ENDDO
228    ENDDO
229
230
231    ! Compute vertical transport
232    DO l=ll_beginp1,ll_end
233       !DIR$ SIMD
234       DO ij=ij_begin,ij_end         
235          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))
236          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))
237          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))
238       ENDDO
239    ENDDO
240
241    !--> flush wwuu 
242    !$OMP BARRIER
243
244    ! Add vertical transport to du
245    DO l=ll_begin,ll_end
246       !DIR$ SIMD
247       DO ij=ij_begin,ij_end         
248          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))
249          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))
250          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))
251       ENDDO
252    ENDDO
253
254    !  DO l=ll_beginp1,ll_end
255    !!DIR$ SIMD
256    !      DO ij=ij_begin,ij_end         
257    !        wwuu_out(ij+u_right,l) = wwuu(ij+u_right,l)
258    !        wwuu_out(ij+u_lup,l)   = wwuu(ij+u_lup,l)
259    !        wwuu_out(ij+u_ldown,l) = wwuu(ij+u_ldown,l)
260    !     ENDDO
261    !   ENDDO
262
263    CALL trace_end("compute_caldyn_vert")
264
265  END SUBROUTINE compute_caldyn_vert
266
267END MODULE caldyn_kernels_base_mod
Note: See TracBrowser for help on using the repository browser.