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

Last change on this file since 403 was 401, checked in by dubos, 8 years ago

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