Context Navigation

caldyn_gcm.f90 @ 519

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

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

Context Navigation

source: codes/icosagcm/trunk/src/caldyn_gcm.f90 @ 519

Download in other formats: