Context Navigation

source: codes/icosagcm/trunk/src/dynamics/caldyn_gcm.F90 @ 580

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

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

Download in other formats: