Context Navigation

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

Last change on this file since 953 was 953, checked in by adurocher, 5 years ago
trunk : GPU implementation with OpenACC ( merge from glcp.idris.fr )
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
[167]	190	USE disvert_mod, ONLY : caldyn_eta, eta_mass
[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)
[387]	321	CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz(:,:,1), du)
[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)
	350	CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
	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.

Download in other formats: