Context Navigation

caldyn_gcm.f90 @ 552

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

Download in other formats: