Context Navigation

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

Last change on this file since 628 was 628, checked in by dubos, 7 years ago
devel : new "conservative" variant for vertical transport of momentum
File size: 12.9 KB

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

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

Download in other formats: