Context Navigation

source: codes/icosagcm/devel/src/dynamics/caldyn_gcm.F90 @ 879

Last change on this file since 879 was 861, checked in by dubos, 5 years ago
devel : compute_X => compute_X_unst or compute_X_hex in init_caldyn
File size: 14.0 KB

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

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

Download in other formats: