Context Navigation

caldyn_gcm.F90 @ 954

Last change on this file since 954 was 954, checked in by adurocher, 5 years ago

trunk : Added metric terms to kernels parameters to avoid Host/GPU transferts

Metric terms are now subroutine parameters instead of module variables in kernel subroutines. Dummy arguments for metric terms are now defined as fixed-size arrays, and arrays dimensions are well known when entering an 'acc data' region. Array descriptors are no longer transferred form host to device each time the data region is executed.

File size: 12.6 KB

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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: codes/icosagcm/trunk/src/dynamics/caldyn_gcm.F90 @ 954

Download in other formats: