Context Navigation

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

Last change on this file since 440 was 413, checked in by dubos, 8 years ago
Improved output
File size: 10.6 KB

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

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

Download in other formats: