source: codes/icosagcm/trunk/src/caldyn_gcm.f90 @ 361

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

ARK2.3 HEVI time stepping - tested with DCMIP4 only
File size: 9.2 KB
Line 
1MODULE caldyn_gcm_mod
2  USE icosa
3  USE transfert_mod
4  USE caldyn_kernels_mod
5  IMPLICIT NONE
6  PRIVATE
7
8  PUBLIC init_caldyn, caldyn_BC, caldyn, req_ps, req_mass
9
10CONTAINS
11 
12  SUBROUTINE init_caldyn
13    USE icosa
14    USE observable_mod
15!    USE exner_mod
16    USE mpipara
17    USE omp_para
18    IMPLICIT NONE
19    CHARACTER(len=255) :: def
20    INTEGER            :: ind
21    REAL(rstd),POINTER :: planetvel(:)
22 
23    def='energy'
24    CALL getin('caldyn_conserv',def)
25    SELECT CASE(TRIM(def))
26    CASE('energy')
27       caldyn_conserv=energy
28    CASE('enstrophy')
29       caldyn_conserv=enstrophy
30    CASE DEFAULT
31       IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_conserv : <', &
32            TRIM(def),'> options are <energy>, <enstrophy>'
33       STOP
34    END SELECT
35    IF (is_master) PRINT *, 'caldyn_conserv=',def
36
37    CALL allocate_caldyn
38
39    DO ind=1,ndomain
40       IF (.NOT. assigned_domain(ind)) CYCLE
41       CALL swap_dimensions(ind)
42       CALL swap_geometry(ind)
43       planetvel=f_planetvel(ind)
44       CALL compute_planetvel(planetvel)
45    END DO
46
47  END SUBROUTINE init_caldyn
48
49  SUBROUTINE allocate_caldyn
50  USE icosa
51  IMPLICIT NONE
52
53    CALL allocate_field(f_out_u,field_u,type_real,llm) 
54    CALL allocate_field(f_qu,field_u,type_real,llm) 
55    CALL allocate_field(f_qv,field_z,type_real,llm) 
56    CALL allocate_field(f_pk,    field_t,type_real,llm,  name='pk')
57    CALL allocate_field(f_wwuu,  field_u,type_real,llm+1,name='wwuu')
58    CALL allocate_field(f_planetvel,  field_u,type_real, name='planetvel') ! planetary velocity at r=a
59
60  END SUBROUTINE allocate_caldyn
61
62  SUBROUTINE caldyn_BC(f_phis, f_geopot, f_wflux)
63    USE icosa
64    USE mpipara
65    USE omp_para
66    TYPE(t_field),POINTER :: f_phis(:)
67    TYPE(t_field),POINTER :: f_geopot(:)
68    TYPE(t_field),POINTER :: f_wflux(:)
69    REAL(rstd),POINTER  :: phis(:)
70    REAL(rstd),POINTER  :: wflux(:,:)
71    REAL(rstd),POINTER  :: geopot(:,:)
72    REAL(rstd),POINTER  :: wwuu(:,:)
73
74    INTEGER :: ind,i,j,ij,l
75
76    IF (is_omp_first_level) THEN
77       DO ind=1,ndomain
78          IF (.NOT. assigned_domain(ind)) CYCLE
79          CALL swap_dimensions(ind)
80          CALL swap_geometry(ind)
81          geopot=f_geopot(ind)
82          phis=f_phis(ind)
83          wflux=f_wflux(ind)
84          wwuu=f_wwuu(ind)
85         
86          DO ij=ij_begin_ext,ij_end_ext
87              ! lower BCs : geopot=phis, wflux=0, wwuu=0
88              geopot(ij,1) = phis(ij)
89              wflux(ij,1) = 0.
90              wwuu(ij+u_right,1)=0   
91              wwuu(ij+u_lup,1)=0   
92              wwuu(ij+u_ldown,1)=0
93              ! top BCs : wflux=0, wwuu=0
94              wflux(ij,llm+1)  = 0.
95              wwuu(ij+u_right,llm+1)=0   
96              wwuu(ij+u_lup,llm+1)=0   
97              wwuu(ij+u_ldown,llm+1)=0
98          ENDDO
99       END DO
100    ENDIF
101
102    !$OMP BARRIER
103  END SUBROUTINE caldyn_BC
104   
105  SUBROUTINE caldyn(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, &
106       f_geopot, f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du)
107    USE icosa
108    USE observable_mod
109    USE disvert_mod, ONLY : caldyn_eta, eta_mass
110    USE vorticity_mod
111    USE kinetic_mod
112    USE theta2theta_rhodz_mod
113    USE wind_mod
114    USE mpipara
115    USE trace
116    USE omp_para
117    USE output_field_mod
118    USE checksum_mod
119    IMPLICIT NONE
120    LOGICAL,INTENT(IN)    :: write_out
121    TYPE(t_field),POINTER :: f_phis(:)
122    TYPE(t_field),POINTER :: f_ps(:)
123    TYPE(t_field),POINTER :: f_mass(:)
124    TYPE(t_field),POINTER :: f_theta_rhodz(:)
125    TYPE(t_field),POINTER :: f_u(:)
126    TYPE(t_field),POINTER :: f_q(:)
127    TYPE(t_field),POINTER :: f_geopot(:)
128    TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:)
129    TYPE(t_field) :: f_dps(:)
130    TYPE(t_field) :: f_dmass(:)
131    TYPE(t_field) :: f_dtheta_rhodz(:)
132    TYPE(t_field) :: f_du(:)
133   
134    REAL(rstd),POINTER :: ps(:), dps(:)
135    REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:), dtheta_rhodz(:,:)
136    REAL(rstd),POINTER :: u(:,:), du(:,:), hflux(:,:), wflux(:,:)
137    REAL(rstd),POINTER :: qu(:,:)
138    REAL(rstd),POINTER :: qv(:,:)
139
140! temporary shared variable
141    REAL(rstd),POINTER  :: theta(:,:) 
142    REAL(rstd),POINTER  :: pk(:,:)
143    REAL(rstd),POINTER  :: geopot(:,:)
144    REAL(rstd),POINTER  :: convm(:,:) 
145    REAL(rstd),POINTER  :: wwuu(:,:)
146       
147    INTEGER :: ind
148    LOGICAL,SAVE :: first=.TRUE.
149!$OMP THREADPRIVATE(first)
150   
151    IF (first) THEN
152      first=.FALSE.
153      IF(caldyn_eta==eta_mass) THEN
154         CALL init_message(f_ps,req_i1,req_ps)
155      ELSE
156         CALL init_message(f_mass,req_i1,req_mass)
157      END IF
158      CALL init_message(f_theta_rhodz,req_i1,req_theta_rhodz)
159      CALL init_message(f_u,req_e1_vect,req_u)
160      CALL init_message(f_qu,req_e1_scal,req_qu)
161      ! Overlapping com/compute (deactivated) : MPI messages need to be sent at first call to caldyn
162      ! This is needed only once : the next ones will be sent by timeloop
163!      IF(caldyn_eta==eta_mass) THEN
164!         CALL send_message(f_ps,req_ps)
165!         CALL wait_message(req_ps) 
166!      ELSE
167!         CALL send_message(f_mass,req_mass)
168!         CALL wait_message(req_mass) 
169!      END IF
170    ENDIF
171   
172    CALL trace_start("caldyn")
173
174    IF(caldyn_eta==eta_mass) THEN
175       CALL send_message(f_ps,req_ps) ! COM00
176    ELSE
177       CALL send_message(f_mass,req_mass) ! COM00
178    END IF
179   
180    CALL send_message(f_theta_rhodz,req_theta_rhodz) ! COM01
181    CALL send_message(f_u,req_u) ! COM02
182
183    SELECT CASE(caldyn_conserv)
184    CASE(energy) ! energy-conserving
185       DO ind=1,ndomain
186          IF (.NOT. assigned_domain(ind)) CYCLE
187          CALL swap_dimensions(ind)
188          CALL swap_geometry(ind)
189          ps=f_ps(ind)
190          u=f_u(ind)
191          theta_rhodz = f_theta_rhodz(ind)
192          mass=f_mass(ind)
193          theta = f_theta(ind)
194          qu=f_qu(ind)
195          qv=f_qv(ind)
196          CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu,qv) ! COM00 COM01 COM02
197       ENDDO
198!       CALL checksum(f_mass)
199!       CALL checksum(f_theta)
200
201       CALL send_message(f_qu,req_qu) ! COM03 wait_message in caldyn_horiz
202!       CALL wait_message(req_qu)
203
204       DO ind=1,ndomain
205          IF (.NOT. assigned_domain(ind)) CYCLE
206          CALL swap_dimensions(ind)
207          CALL swap_geometry(ind)
208          ps=f_ps(ind)
209          u=f_u(ind)
210          theta_rhodz=f_theta_rhodz(ind)
211          mass=f_mass(ind)
212          theta = f_theta(ind)
213          qu=f_qu(ind)
214          pk = f_pk(ind)
215          geopot = f_geopot(ind) 
216          CALL compute_geopot(ps,mass,theta, pk,geopot)
217          hflux=f_hflux(ind)
218          convm = f_dmass(ind)
219          dtheta_rhodz=f_dtheta_rhodz(ind)
220          du=f_du(ind)
221          CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz,du)
222          IF(caldyn_eta==eta_mass) THEN
223             wflux=f_wflux(ind)
224             wwuu=f_wwuu(ind)
225             dps=f_dps(ind)
226             CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
227          END IF
228       ENDDO       
229   
230!       CALL checksum(f_geopot)
231!       CALL checksum(f_dmass)
232!       CALL checksum(f_pk)
233!       CALL checksum(f_pk)
234       
235    CASE(enstrophy) ! enstrophy-conserving
236       DO ind=1,ndomain
237          IF (.NOT. assigned_domain(ind)) CYCLE
238          CALL swap_dimensions(ind)
239          CALL swap_geometry(ind)
240          ps=f_ps(ind)
241          u=f_u(ind)
242          theta_rhodz=f_theta_rhodz(ind)
243          mass=f_mass(ind)
244          theta = f_theta(ind)
245          qu=f_qu(ind)
246          qv=f_qv(ind)
247          CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu,qv)
248          pk = f_pk(ind)
249          geopot = f_geopot(ind) 
250          CALL compute_geopot(ps,mass,theta, pk,geopot)
251          hflux=f_hflux(ind)
252          convm = f_dmass(ind)
253          dtheta_rhodz=f_dtheta_rhodz(ind)
254          du=f_du(ind)
255          CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz,du)
256          IF(caldyn_eta==eta_mass) THEN
257             wflux=f_wflux(ind)
258             wwuu=f_wwuu(ind)
259             dps=f_dps(ind)
260             CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
261          END IF
262       ENDDO
263       
264    CASE DEFAULT
265       STOP
266    END SELECT
267
268!$OMP BARRIER
269    !    CALL check_mass_conservation(f_ps,f_dps)
270    CALL trace_end("caldyn")
271!!$OMP BARRIER
272   
273END SUBROUTINE caldyn
274
275!-------------------------------- Diagnostics ----------------------------
276
277  SUBROUTINE check_mass_conservation(f_ps,f_dps)
278  USE icosa
279  USE mpipara
280  IMPLICIT NONE
281    TYPE(t_field),POINTER :: f_ps(:)
282    TYPE(t_field),POINTER :: f_dps(:)
283    REAL(rstd),POINTER :: ps(:)
284    REAL(rstd),POINTER :: dps(:)
285    REAL(rstd) :: mass_tot,dmass_tot
286    INTEGER :: ind,i,j,ij
287   
288    mass_tot=0
289    dmass_tot=0
290   
291    CALL transfert_request(f_dps,req_i1)
292    CALL transfert_request(f_ps,req_i1)
293
294    DO ind=1,ndomain
295      CALL swap_dimensions(ind)
296      CALL swap_geometry(ind)
297
298      ps=f_ps(ind)
299      dps=f_dps(ind)
300
301      DO j=jj_begin,jj_end
302        DO i=ii_begin,ii_end
303          ij=(j-1)*iim+i
304          IF (domain(ind)%own(i,j)) THEN
305            mass_tot=mass_tot+ps(ij)*Ai(ij)/g
306            dmass_tot=dmass_tot+dps(ij)*Ai(ij)/g
307          ENDIF
308        ENDDO
309      ENDDO
310   
311    ENDDO
312    IF (is_mpi_root) PRINT*, "mass_tot ", mass_tot,"      dmass_tot ",dmass_tot       
313
314  END SUBROUTINE check_mass_conservation 
315 
316END MODULE caldyn_gcm_mod
Note: See TracBrowser for help on using the repository browser.