source: codes/icosagcm/trunk/src/dynamics/caldyn_gcm.f90 @ 552

Last change on this file since 552 was 552, checked in by dubos, 7 years ago

trunk : backported commit 534 from devel
File size: 11.8 KB
Line 
1MODULE caldyn_gcm_mod
2  USE icosa
3  USE transfert_mod
4  USE caldyn_kernels_hevi_mod
5  USE caldyn_kernels_base_mod
6  USE caldyn_kernels_mod
7  IMPLICIT NONE
8  PRIVATE
9
10  PUBLIC init_caldyn, caldyn_BC, caldyn
11 
12CONTAINS
13 
14  SUBROUTINE init_caldyn
15    USE icosa
16    USE observable_mod
17    USE mpipara
18    USE omp_para
19    IMPLICIT NONE
20    CHARACTER(len=255) :: def
21    INTEGER            :: ind
22    REAL(rstd),POINTER :: planetvel(:)
23
24    hydrostatic=.TRUE.
25    CALL getin("hydrostatic",hydrostatic)
26 
27    def='energy'
28    CALL getin('caldyn_conserv',def)
29    SELECT CASE(TRIM(def))
30    CASE('energy')
31       caldyn_conserv=energy
32    CASE('enstrophy')
33       caldyn_conserv=enstrophy
34    CASE DEFAULT
35       IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_conserv : <', &
36            TRIM(def),'> options are <energy>, <enstrophy>'
37       STOP
38    END SELECT
39    IF (is_master) PRINT *, 'caldyn_conserv=',def
40
41    nqdyn=1 ! default value
42    physics_thermo = thermo_none
43
44    def='theta'
45    CALL getin('thermo',def)
46    SELECT CASE(TRIM(def))
47    CASE('boussinesq')
48       boussinesq=.TRUE.
49       caldyn_thermo=thermo_boussinesq
50       IF(.NOT. hydrostatic) THEN
51          PRINT *, 'thermo=boussinesq and hydrostatic=.FALSE. : Non-hydrostatic boussinesq equations are not supported'
52          STOP
53       END IF
54    CASE('theta')
55       caldyn_thermo=thermo_theta
56       physics_thermo=thermo_dry
57    CASE('entropy')
58       caldyn_thermo=thermo_entropy
59       physics_thermo=thermo_dry
60    CASE('theta_fake_moist')
61       caldyn_thermo=thermo_theta
62       physics_thermo=thermo_fake_moist
63    CASE('entropy_fake_moist')
64       caldyn_thermo=thermo_entropy
65       physics_thermo=thermo_fake_moist
66    CASE('moist')
67       caldyn_thermo=thermo_moist_debug
68       physics_thermo=thermo_moist
69       nqdyn = 2
70    CASE DEFAULT
71       IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_thermo : <', &
72            TRIM(def),'> options are <theta>, <entropy>'
73       STOP
74    END SELECT
75
76    IF(is_master) THEN
77       SELECT CASE(caldyn_thermo)
78       CASE(thermo_theta)
79          PRINT *, 'caldyn_thermo = thermo_theta'
80       CASE(thermo_entropy)
81          PRINT *, 'caldyn_thermo = thermo_entropy'
82       CASE(thermo_moist_debug)
83          PRINT *, 'caldyn_thermo = thermo_moist_debug'
84       CASE DEFAULT
85          STOP
86       END SELECT
87
88       SELECT CASE(physics_thermo)
89       CASE(thermo_dry)
90          PRINT *, 'physics_thermo = thermo_dry'
91       CASE(thermo_fake_moist)
92          PRINT *, 'physics_thermo = thermo_fake_moist'
93       CASE(thermo_moist)
94          PRINT *, 'physics_thermo = thermo_moist'
95       END SELECT
96
97       PRINT *, 'nqdyn =', nqdyn
98    END IF
99
100    CALL allocate_caldyn
101
102    DO ind=1,ndomain
103       IF (.NOT. assigned_domain(ind)) CYCLE
104       CALL swap_dimensions(ind)
105       CALL swap_geometry(ind)
106       planetvel=f_planetvel(ind)
107       CALL compute_planetvel(planetvel)
108    END DO
109
110  END SUBROUTINE init_caldyn
111
112  SUBROUTINE allocate_caldyn
113  USE icosa
114  IMPLICIT NONE
115
116    CALL allocate_field(f_out_u,field_u,type_real,llm) 
117    CALL allocate_field(f_qu,field_u,type_real,llm) 
118    CALL allocate_field(f_qv,field_z,type_real,llm) 
119    CALL allocate_field(f_pk,    field_t,type_real,llm,  name='pk')
120    CALL allocate_field(f_wwuu,  field_u,type_real,llm+1,name='wwuu')
121    CALL allocate_field(f_planetvel,  field_u,type_real, name='planetvel') ! planetary velocity at r=a
122
123  END SUBROUTINE allocate_caldyn
124
125  SUBROUTINE caldyn_BC(f_phis, f_geopot, f_wflux)
126    USE icosa
127    USE mpipara
128    USE omp_para
129    TYPE(t_field),POINTER :: f_phis(:)
130    TYPE(t_field),POINTER :: f_geopot(:)
131    TYPE(t_field),POINTER :: f_wflux(:)
132    REAL(rstd),POINTER  :: phis(:)
133    REAL(rstd),POINTER  :: wflux(:,:)
134    REAL(rstd),POINTER  :: geopot(:,:)
135    REAL(rstd),POINTER  :: wwuu(:,:)
136
137    INTEGER :: ind,i,j,ij,l
138
139    IF (is_omp_first_level) THEN
140       DO ind=1,ndomain
141          IF (.NOT. assigned_domain(ind)) CYCLE
142          CALL swap_dimensions(ind)
143          CALL swap_geometry(ind)
144          geopot=f_geopot(ind)
145          phis=f_phis(ind)
146          wflux=f_wflux(ind)
147          wwuu=f_wwuu(ind)
148         
149          DO ij=ij_begin_ext,ij_end_ext
150              ! lower BCs : geopot=phis, wflux=0, wwuu=0
151              geopot(ij,1) = phis(ij)
152              wflux(ij,1) = 0.
153              wwuu(ij+u_right,1)=0   
154              wwuu(ij+u_lup,1)=0   
155              wwuu(ij+u_ldown,1)=0
156              ! top BCs : wflux=0, wwuu=0
157              wflux(ij,llm+1)  = 0.
158              wwuu(ij+u_right,llm+1)=0   
159              wwuu(ij+u_lup,llm+1)=0   
160              wwuu(ij+u_ldown,llm+1)=0
161          ENDDO
162       END DO
163    ENDIF
164
165    !$OMP BARRIER
166  END SUBROUTINE caldyn_BC
167   
168  SUBROUTINE caldyn(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, &
169       f_geopot, f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du)
170    USE icosa
171    USE observable_mod
172    USE disvert_mod, ONLY : caldyn_eta, eta_mass
173    USE vorticity_mod
174    USE kinetic_mod
175    USE theta2theta_rhodz_mod
176    USE wind_mod
177    USE mpipara
178    USE trace
179    USE omp_para
180    USE output_field_mod
181    USE checksum_mod
182    IMPLICIT NONE
183    LOGICAL,INTENT(IN)    :: write_out
184    TYPE(t_field),POINTER :: f_phis(:)
185    TYPE(t_field),POINTER :: f_ps(:)
186    TYPE(t_field),POINTER :: f_mass(:)
187    TYPE(t_field),POINTER :: f_theta_rhodz(:)
188    TYPE(t_field),POINTER :: f_u(:)
189    TYPE(t_field),POINTER :: f_q(:)
190    TYPE(t_field),POINTER :: f_geopot(:)
191    TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:)
192    TYPE(t_field) :: f_dps(:)
193    TYPE(t_field) :: f_dmass(:)
194    TYPE(t_field) :: f_dtheta_rhodz(:)
195    TYPE(t_field) :: f_du(:)
196   
197    REAL(rstd),POINTER :: ps(:), dps(:)
198    REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:,:), dtheta_rhodz(:,:,:)
199    REAL(rstd),POINTER :: u(:,:), du(:,:), hflux(:,:), wflux(:,:)
200    REAL(rstd),POINTER :: qu(:,:)
201    REAL(rstd),POINTER :: qv(:,:)
202
203! temporary shared variable
204    REAL(rstd),POINTER  :: theta(:,:,:) 
205    REAL(rstd),POINTER  :: pk(:,:)
206    REAL(rstd),POINTER  :: geopot(:,:)
207    REAL(rstd),POINTER  :: convm(:,:) 
208    REAL(rstd),POINTER  :: wwuu(:,:)
209       
210    INTEGER :: ind
211    LOGICAL,SAVE :: first=.TRUE.
212!$OMP THREADPRIVATE(first)
213   
214    IF (first) THEN
215      first=.FALSE.
216      IF(caldyn_eta==eta_mass) THEN
217         CALL init_message(f_ps,req_i1,req_ps)
218      ELSE
219         CALL init_message(f_mass,req_i1,req_mass)
220      END IF
221      CALL init_message(f_theta_rhodz,req_i1,req_theta_rhodz)
222      CALL init_message(f_u,req_e1_vect,req_u)
223      CALL init_message(f_qu,req_e1_scal,req_qu)
224      ! Overlapping com/compute (deactivated) : MPI messages need to be sent at first call to caldyn
225      ! This is needed only once : the next ones will be sent by timeloop
226!      IF(caldyn_eta==eta_mass) THEN
227!         CALL send_message(f_ps,req_ps)
228!         CALL wait_message(req_ps) 
229!      ELSE
230!         CALL send_message(f_mass,req_mass)
231!         CALL wait_message(req_mass) 
232!      END IF
233    ENDIF
234   
235    CALL trace_start("caldyn")
236
237    IF(caldyn_eta==eta_mass) THEN
238       CALL send_message(f_ps,req_ps) ! COM00
239       CALL wait_message(req_ps) ! COM00
240    ELSE
241       CALL send_message(f_mass,req_mass) ! COM00
242       CALL wait_message(req_mass) ! COM00
243    END IF
244   
245    CALL send_message(f_theta_rhodz,req_theta_rhodz) ! COM01
246    CALL wait_message(req_theta_rhodz) ! COM01 Moved from caldyn_pvort
247    CALL send_message(f_u,req_u) ! COM02
248    CALL wait_message(req_u) ! COM02
249
250    IF(.NOT.hydrostatic) THEN
251       STOP 'caldyn_gcm may not be used yet when non-hydrostatic'
252    END IF
253
254    SELECT CASE(caldyn_conserv)
255    CASE(energy) ! energy-conserving
256       DO ind=1,ndomain
257          IF (.NOT. assigned_domain(ind)) CYCLE
258          CALL swap_dimensions(ind)
259          CALL swap_geometry(ind)
260          ps=f_ps(ind)
261          u=f_u(ind)
262          theta_rhodz = f_theta_rhodz(ind)
263          mass=f_mass(ind)
264          theta = f_theta(ind)
265          qu=f_qu(ind)
266          qv=f_qv(ind)
267          pk = f_pk(ind)
268          geopot = f_geopot(ind) 
269          hflux=f_hflux(ind)
270          convm = f_dmass(ind)
271          dtheta_rhodz=f_dtheta_rhodz(ind)
272          du=f_du(ind)
273          CALL compute_pvort(ps,u,theta_rhodz(:,:,1), mass,theta,qu,qv) ! COM00 COM01 COM02
274!          CALL compute_theta(ps,theta_rhodz, mass,theta)
275!          CALL compute_pvort_only(u,mass,qu,qv)
276
277          CALL compute_geopot(mass,theta, ps,pk,geopot)
278!          du(:,:)=0.
279!          CALL compute_caldyn_fast(0.,u,mass,theta,pk,geopot,du)
280       ENDDO       
281
282       CALL send_message(f_u,req_u) ! COM02
283       CALL wait_message(req_u) ! COM02
284       CALL send_message(f_qu,req_qu) ! COM03
285       CALL wait_message(req_qu) ! COM03
286
287       DO ind=1,ndomain
288          IF (.NOT. assigned_domain(ind)) CYCLE
289          CALL swap_dimensions(ind)
290          CALL swap_geometry(ind)
291          ps=f_ps(ind)
292          u=f_u(ind)
293          theta_rhodz = f_theta_rhodz(ind)
294          mass=f_mass(ind)
295          theta = f_theta(ind)
296          qu=f_qu(ind)
297          qv=f_qv(ind)
298          pk = f_pk(ind)
299          geopot = f_geopot(ind) 
300          hflux=f_hflux(ind)
301          convm = f_dmass(ind)
302          dtheta_rhodz=f_dtheta_rhodz(ind)
303          du=f_du(ind)
304
305          CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz(:,:,1),du)
306!          CALL compute_caldyn_slow_hydro(u,mass,hflux,du, .FALSE.) ! FIXME
307!          CALL compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du)
308          IF(caldyn_eta==eta_mass) THEN
309             wflux=f_wflux(ind)
310             wwuu=f_wwuu(ind)
311             dps=f_dps(ind)
312             CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz(:,:,1), du)
313          END IF
314       ENDDO       
315   
316    CASE(enstrophy) ! enstrophy-conserving
317       DO ind=1,ndomain
318          IF (.NOT. assigned_domain(ind)) CYCLE
319          CALL swap_dimensions(ind)
320          CALL swap_geometry(ind)
321          ps=f_ps(ind)
322          u=f_u(ind)
323          theta_rhodz=f_theta_rhodz(ind)
324          mass=f_mass(ind)
325          theta = f_theta(ind)
326          qu=f_qu(ind)
327          qv=f_qv(ind)
328          CALL compute_pvort(ps,u,theta_rhodz(:,:,1), mass,theta,qu,qv)
329          pk = f_pk(ind)
330          geopot = f_geopot(ind) 
331          CALL compute_geopot(ps,mass,theta, pk,geopot)
332          hflux=f_hflux(ind)
333          convm = f_dmass(ind)
334          dtheta_rhodz=f_dtheta_rhodz(ind)
335          du=f_du(ind)
336          CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz(:,:,1),du)
337          IF(caldyn_eta==eta_mass) THEN
338             wflux=f_wflux(ind)
339             wwuu=f_wwuu(ind)
340             dps=f_dps(ind)
341             CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
342          END IF
343       ENDDO
344       
345    CASE DEFAULT
346       STOP
347    END SELECT
348
349!$OMP BARRIER
350    !    CALL check_mass_conservation(f_ps,f_dps)
351    CALL trace_end("caldyn")
352!!$OMP BARRIER
353   
354END SUBROUTINE caldyn
355
356!-------------------------------- Diagnostics ----------------------------
357
358  SUBROUTINE check_mass_conservation(f_ps,f_dps)
359  USE icosa
360  USE mpipara
361  IMPLICIT NONE
362    TYPE(t_field),POINTER :: f_ps(:)
363    TYPE(t_field),POINTER :: f_dps(:)
364    REAL(rstd),POINTER :: ps(:)
365    REAL(rstd),POINTER :: dps(:)
366    REAL(rstd) :: mass_tot,dmass_tot
367    INTEGER :: ind,i,j,ij
368   
369    mass_tot=0
370    dmass_tot=0
371   
372    CALL transfert_request(f_dps,req_i1)
373    CALL transfert_request(f_ps,req_i1)
374
375    DO ind=1,ndomain
376      CALL swap_dimensions(ind)
377      CALL swap_geometry(ind)
378
379      ps=f_ps(ind)
380      dps=f_dps(ind)
381
382      DO j=jj_begin,jj_end
383        DO i=ii_begin,ii_end
384          ij=(j-1)*iim+i
385          IF (domain(ind)%own(i,j)) THEN
386            mass_tot=mass_tot+ps(ij)*Ai(ij)/g
387            dmass_tot=dmass_tot+dps(ij)*Ai(ij)/g
388          ENDIF
389        ENDDO
390      ENDDO
391   
392    ENDDO
393    IF (is_mpi_root) PRINT*, "mass_tot ", mass_tot,"      dmass_tot ",dmass_tot       
394
395  END SUBROUTINE check_mass_conservation 
396 
397END MODULE caldyn_gcm_mod
Note: See TracBrowser for help on using the repository browser.