Context Navigation

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

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

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

Download in other formats: