Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

physics.f90

Last change on this file was 1029, checked in by dubos, 4 years ago
devel : fix iteration number passed to physics_column
File size: 13.6 KB

Rev	Line
[81]	1	MODULE physics_mod
[714]	2	USE icosa
[196]	3	USE field_mod
[714]	4	USE physics_interface_mod
	5	USE omp_para
	6	IMPLICIT NONE
[196]	7	PRIVATE
	8
[739]	9	INTEGER, PARAMETER :: phys_none=0, phys_column=1, &
	10	phys_HS94=3, phys_LB2012=4, &
[1007]	11	phys_DCMIP=11, phys_DCMIP2016=12, phys_plugin=13, &
[739]	12	phys_lmdz_generic=21, phys_external=22
[213]	13	INTEGER :: phys_type
[481]	14	TYPE(t_field),POINTER,SAVE :: f_dulon(:), f_dulat(:)
	15	TYPE(t_field),POINTER,SAVE :: f_ulon(:), f_ulat(:)
	16	TYPE(t_field),POINTER,SAVE :: f_p(:), f_pk(:)
	17	TYPE(t_field),POINTER,SAVE :: f_temp(:)
[714]	18	TYPE(t_field),POINTER,SAVE :: f_du_phys(:)
[196]	19
[481]	20	CHARACTER(LEN=255),SAVE :: physics_type
[186]	21	!$OMP THREADPRIVATE(physics_type)
[81]	22
[1007]	23	INTERFACE
	24	SUBROUTINE plugin_init_physics
	25	END SUBROUTINE plugin_init_physics
[213]	26
[1007]	27	SUBROUTINE plugin_physics
	28	END SUBROUTINE plugin_physics
	29	END INTERFACE
	30
	31	PROCEDURE(plugin_init_physics), POINTER :: init_physics_plugin => NULL()
	32	PROCEDURE(plugin_physics), POINTER :: physics_plugin => NULL()
	33
	34	PUBLIC :: physics, init_physics, init_physics_plugin, physics_plugin, zero_du_phys
	35
[81]	36	CONTAINS
	37
[98]	38	SUBROUTINE init_physics
[178]	39	USE mpipara
	40	USE etat0_mod
[739]	41	USE etat0_venus_mod, ONLY : init_phys_venus=>init_physics
[215]	42	USE physics_dcmip_mod, ONLY : init_physics_dcmip=>init_physics
[381]	43	USE physics_dcmip2016_mod, ONLY : init_physics_dcmip2016=>init_physics
[327]	44	USE physics_lmdz_generic_mod, ONLY : init_physics_lmdz_generic=>init_physics
[347]	45	USE physics_external_mod, ONLY : init_physics_external=>init_physics
[739]	46	LOGICAL :: done
	47	physics_inout%dt_phys = dt*itau_physics
	48	!$OMP PARALLEL
	49	CALL allocate_field(f_du_phys,field_u,type_real,llm, name='du_phys')
[170]	50
[217]	51	physics_type='none'
[81]	52	CALL getin("physics",physics_type)
[739]	53	! below, flag done is set to .FALSE. if the CALL to init_XXX must be done outside any OMP PARALLEL region
	54	done=.TRUE.
	55	phys_type=phys_column
[81]	56	SELECT CASE(TRIM(physics_type))
[217]	57	CASE ('none')
	58	IF(is_mpi_root) PRINT*,"NO PHYSICAL PACKAGE USED"
	59	phys_type = phys_none
	60	CASE ('held_suarez')
	61	phys_type = phys_HS94
[327]	62	CASE ('phys_lmdz_generic')
	63	phys_type=phys_lmdz_generic
[739]	64	done = .FALSE.
[347]	65	CASE ('phys_external')
	66	phys_type=phys_external
[739]	67	done = .FALSE.
	68	END SELECT
[726]	69
[739]	70	IF(phys_type == phys_column) THEN
[217]	71	CALL allocate_field(f_dulon,field_t,type_real,llm, name='dulon')
	72	CALL allocate_field(f_dulat,field_t,type_real,llm, name='dulat')
[295]	73	CALL allocate_field(f_temp,field_t,type_real,llm, name='temp')
[481]	74	CALL allocate_field(f_ulon,field_t,type_real,llm, name='ulon')
	75	CALL allocate_field(f_ulat,field_t,type_real,llm, name='ulat')
	76	CALL allocate_field(f_p,field_t,type_real,llm+1, name='p')
	77	CALL allocate_field(f_pk,field_t,type_real,llm, name='pk')
[217]	78	CALL init_pack_before ! Compute physics_inout%ngrid and offsets used by pack/unpack
[739]	79	CALL init_pack_after ! Defines Ai, lon, lat in physics_inout
	80
	81	SELECT CASE(TRIM(physics_type))
	82	CASE ('dcmip')
	83	phys_type = phys_DCMIP
	84	CALL init_physics_dcmip
	85	CASE ('dcmip2016')
	86	phys_type = phys_DCMIP2016
	87	CALL init_physics_dcmip2016
[741]	88	CASE ('Lebonnois2012')
	89	phys_type = phys_LB2012
	90	CALL init_phys_venus
[1007]	91	CASE ('plugin')
	92	IF(.NOT.ASSOCIATED(init_physics_plugin)) THEN
	93	PRINT *, 'FATAL : physics = plugin selected by init_physics_plugin not set by driver program. Aborting.'
	94	STOP
	95	END IF
	96	phys_type = phys_plugin
	97	CALL init_physics_plugin
[739]	98	CASE DEFAULT
	99	IF(is_mpi_root) PRINT*, 'init_physics : Bad selector for variable physics <',&
	100	TRIM(physics_type), '> options are <none>, <held_suarez>, <Lebonnois2012>,', &
[1007]	101	'<dcmip>, <dcmip2016>, <plugin>, <phys_lmdz_generic>, <phys_external>'
[739]	102	STOP
	103	END SELECT
	104
	105	END IF
[726]	106	!$OMP END PARALLEL
	107
[1026]	108	IF(.NOT. done) THEN
[739]	109	SELECT CASE(phys_type)
	110	CASE(phys_external)
	111	CALL init_physics_external
	112	CASE(phys_lmdz_generic)
	113	CALL init_physics_lmdz_generic
	114	END SELECT
	115	END IF
[726]	116
[739]	117	IF(is_mpi_root) PRINT *, 'phys_type = ',phys_type
[726]	118
[81]	119	END SUBROUTINE init_physics
	120
[714]	121	SUBROUTINE zero_du_phys()
	122	REAL(rstd), DIMENSION(:,:), POINTER :: du
	123	INTEGER :: ind
	124	DO ind=1,ndomain
	125	IF (.NOT. assigned_domain(ind)) CYCLE
	126	CALL swap_dimensions(ind)
	127	CALL swap_geometry(ind)
	128	du=f_du_phys(ind)
	129	du(:,ll_begin:ll_end) = 0.
	130	END DO
	131	END SUBROUTINE zero_du_phys
	132
	133	SUBROUTINE add_du_phys(coef, f_u)
	134	REAL(rstd), INTENT(IN) :: coef ! -1 before physics, +1 after physics
	135	TYPE(t_field),POINTER :: f_u(:) ! velocity field before/after call to physics
	136	REAL(rstd), DIMENSION(:,:), POINTER :: u, du
	137	INTEGER :: ind
	138	DO ind=1,ndomain
	139	IF (.NOT. assigned_domain(ind)) CYCLE
	140	CALL swap_dimensions(ind)
	141	CALL swap_geometry(ind)
	142	du=f_du_phys(ind)
	143	u=f_u(ind)
	144	du(:,ll_begin:ll_end) = du(:,ll_begin:ll_end) + coef*u(:,ll_begin:ll_end)
	145	END DO
	146	END SUBROUTINE add_du_phys
	147
[741]	148	SUBROUTINE physics(it,f_phis, f_geopot, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q)
[327]	149	USE physics_lmdz_generic_mod, ONLY : physics_lmdz_generic => physics
[347]	150	USE physics_external_mod, ONLY : physics_external => physics
[213]	151	USE physics_dcmip_mod, ONLY : write_physics_dcmip => write_physics
[381]	152	USE physics_dcmip2016_mod, ONLY : write_physics_dcmip2016 => write_physics
[170]	153	USE etat0_heldsz_mod
[99]	154	INTEGER, INTENT(IN) :: it
[81]	155	TYPE(t_field),POINTER :: f_phis(:)
[741]	156	TYPE(t_field),POINTER :: f_geopot(:)
[81]	157	TYPE(t_field),POINTER :: f_ps(:)
	158	TYPE(t_field),POINTER :: f_theta_rhodz(:)
	159	TYPE(t_field),POINTER :: f_ue(:)
[327]	160	TYPE(t_field),POINTER :: f_wflux(:)
[81]	161	TYPE(t_field),POINTER :: f_q(:)
[196]	162
[149]	163	LOGICAL:: firstcall,lastcall
[196]	164	INTEGER :: ind
[213]	165	TYPE(t_physics_inout) :: args
[170]	166
[714]	167	IF(MOD(it,itau_physics)==0 .AND. phys_type/=phys_none) THEN
	168
	169	! as a result of the the two calls to add_du_phys,
	170	! du_phys increases by u(after physics) - u (before physics)
	171	CALL add_du_phys(-1., f_ue)
[1011]	172
[213]	173	SELECT CASE(phys_type)
	174	CASE(phys_HS94)
	175	CALL held_suarez(f_ps,f_theta_rhodz,f_ue)
[327]	176	CASE (phys_lmdz_generic)
	177	CALL physics_lmdz_generic(it ,f_phis, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q)
[347]	178	CASE (phys_external)
	179	CALL physics_external(it ,f_phis, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q)
[213]	180	CASE DEFAULT
[1011]	181	CALL transfert_request(f_ue, req_e1_vect)
[1029]	182	CALL physics_column(it/itau_physics, f_phis, f_geopot, f_ps, f_theta_rhodz, f_ue, f_q)
[213]	183	END SELECT
[149]	184
[213]	185	CALL transfert_request(f_theta_rhodz,req_i0)
	186	CALL transfert_request(f_ue,req_e0_vect)
	187	CALL transfert_request(f_q,req_i0)
[714]	188
	189	CALL add_du_phys(1., f_ue)
[213]	190	END IF
[149]	191
[213]	192	IF (mod(it,itau_out)==0 ) THEN
[714]	193	CALL write_physics_tendencies
	194	CALL zero_du_phys
[213]	195	SELECT CASE(phys_type)
	196	CASE (phys_DCMIP)
	197	CALL write_physics_dcmip
[381]	198	CASE (phys_DCMIP2016)
	199	CALL write_physics_dcmip2016
[213]	200	END SELECT
	201	END IF
	202
	203	END SUBROUTINE physics
	204
[714]	205	SUBROUTINE write_physics_tendencies
	206	USE observable_mod, ONLY : f_buf_ulon, f_buf_ulat
	207	USE wind_mod
	208	USE output_field_mod
	209	CALL transfert_request(f_du_phys,req_e1_vect)
	210	CALL un2ulonlat(f_du_phys, f_buf_ulon, f_buf_ulat, (1./(dt*itau_out)))
	211	CALL output_field("dulon_phys",f_buf_ulon)
	212	CALL output_field("dulat_phys",f_buf_ulat)
	213	END SUBROUTINE write_physics_tendencies
	214
[741]	215	SUBROUTINE physics_column(it, f_phis, f_geopot, f_ps, f_theta_rhodz, f_ue, f_q)
[213]	216	USE physics_dcmip_mod, ONLY : full_physics_dcmip => full_physics
[381]	217	USE physics_dcmip2016_mod, ONLY : full_physics_dcmip2016 => full_physics
[741]	218	USE etat0_venus_mod, ONLY : full_physics_venus=>full_physics
[295]	219	USE theta2theta_rhodz_mod
[281]	220	USE mpipara
[481]	221	USE checksum_mod
[213]	222	TYPE(t_field),POINTER :: f_phis(:)
[741]	223	TYPE(t_field),POINTER :: f_geopot(:)
[213]	224	TYPE(t_field),POINTER :: f_ps(:)
	225	TYPE(t_field),POINTER :: f_theta_rhodz(:)
	226	TYPE(t_field),POINTER :: f_ue(:)
	227	TYPE(t_field),POINTER :: f_q(:)
	228	REAL(rstd),POINTER :: phis(:)
[741]	229	REAL(rstd),POINTER :: geopot(:,:)
[213]	230	REAL(rstd),POINTER :: ps(:)
[295]	231	REAL(rstd),POINTER :: temp(:,:)
[213]	232	REAL(rstd),POINTER :: ue(:,:)
	233	REAL(rstd),POINTER :: dulon(:,:)
	234	REAL(rstd),POINTER :: dulat(:,:)
	235	REAL(rstd),POINTER :: q(:,:,:)
[481]	236	REAL(rstd),POINTER :: p(:,:)
	237	REAL(rstd),POINTER :: pk(:,:)
	238	REAL(rstd),POINTER :: ulon(:,:)
	239	REAL(rstd),POINTER :: ulat(:,:)
[213]	240	INTEGER :: it, ind
	241
[1011]	242	physics_inout%it = it
	243
[295]	244	CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_temp)
[381]	245
[213]	246	DO ind=1,ndomain
	247	IF (.NOT. assigned_domain(ind)) CYCLE
	248	CALL swap_dimensions(ind)
	249	CALL swap_geometry(ind)
	250	phis=f_phis(ind)
[741]	251	geopot=f_geopot(ind)
[213]	252	ps=f_ps(ind)
[295]	253	temp=f_temp(ind)
[213]	254	ue=f_ue(ind)
	255	q=f_q(ind)
[481]	256	p=f_p(ind)
	257	pk=f_pk(ind)
	258	ulon=f_ulon(ind)
	259	ulat=f_ulat(ind)
[741]	260	CALL pack_physics(pack_info(ind), phis, geopot, ps, temp, ue, q, p, pk, ulon, ulat)
[213]	261	END DO
	262
	263	SELECT CASE(phys_type)
	264	CASE (phys_DCMIP)
[481]	265	IF (is_omp_level_master) CALL full_physics_dcmip
[381]	266	CASE (phys_DCMIP2016)
[481]	267	IF (is_omp_level_master) CALL full_physics_dcmip2016
[741]	268	CASE(phys_LB2012)
	269	IF (is_omp_level_master) CALL full_physics_venus
[1007]	270	CASE(phys_plugin)
	271	IF (is_omp_level_master) CALL physics_plugin
[281]	272	CASE DEFAULT
[481]	273	IF(is_master) PRINT *,'Internal error : illegal value of phys_type', phys_type
[281]	274	STOP
[196]	275	END SELECT
[170]	276
[213]	277	DO ind=1,ndomain
	278	IF (.NOT. assigned_domain(ind)) CYCLE
	279	CALL swap_dimensions(ind)
	280	CALL swap_geometry(ind)
	281	ps=f_ps(ind)
[295]	282	temp=f_temp(ind)
[213]	283	q=f_q(ind)
	284	dulon=f_dulon(ind)
	285	dulat=f_dulat(ind)
[1011]	286	! unpack_physics updates temp and q ; a transfer is required before updating u
[295]	287	CALL unpack_physics(pack_info(ind), ps, temp, q, dulon, dulat)
[213]	288	END DO
[481]	289
[295]	290	CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz)
[170]	291
[213]	292	! Transfer dulon, dulat
	293	CALL transfert_request(f_dulon,req_i0)
	294	CALL transfert_request(f_dulat,req_i0)
	295
	296	DO ind=1,ndomain
	297	IF (.NOT. assigned_domain(ind)) CYCLE
	298	CALL swap_dimensions(ind)
	299	CALL swap_geometry(ind)
	300	ue=f_ue(ind)
	301	dulon=f_dulon(ind)
	302	dulat=f_dulat(ind)
	303	CALL compute_update_velocity(dulon, dulat, ue)
	304	END DO
	305
[214]	306	END SUBROUTINE physics_column
[213]	307
[741]	308	SUBROUTINE pack_physics(info, phis, geopot, ps, temp, ue, q, p, pk, ulon, ulat )
[213]	309	USE wind_mod
[913]	310	USE compute_diagnostics_mod
[213]	311	USE theta2theta_rhodz_mod
[381]	312	USE exner_mod
[213]	313	TYPE(t_pack_info) :: info
	314	REAL(rstd) :: phis(iim*jjm)
[741]	315	REAL(rstd) :: geopot(iim*jjm,llm+1)
[213]	316	REAL(rstd) :: ps(iim*jjm)
[295]	317	REAL(rstd) :: temp(iim*jjm,llm)
[381]	318	REAL(rstd) :: pks(iim*jjm)
	319	REAL(rstd) :: pk(iim*jjm,llm)
[213]	320	REAL(rstd) :: ue(3iimjjm,llm)
	321	REAL(rstd) :: q(iim*jjm,llm,nqtot)
	322
	323	REAL(rstd) :: p(iim*jjm,llm+1)
[584]	324	REAL(rstd) :: uc(iim*jjm,llm,3)
[213]	325	REAL(rstd) :: ulon(iim*jjm,llm)
	326	REAL(rstd) :: ulat(iim*jjm,llm)
	327
[295]	328	!$OMP BARRIER
[1011]	329	CALL compute_pression(ps,p,0) ! pressure at interfaces
[295]	330	!$OMP BARRIER
[381]	331	CALL compute_exner(ps,p,pks,pk,0)
	332	!$OMP BARRIER
[213]	333	CALL compute_wind_centered(ue,uc)
	334	CALL compute_wind_centered_lonlat_compound(uc, ulon, ulat)
[481]	335	!$OMP BARRIER
	336	IF (is_omp_level_master) THEN
	337	CALL pack_domain(info, phis, physics_inout%phis)
[741]	338	CALL pack_domain(info, geopot, physics_inout%geopot)
[481]	339	CALL pack_domain(info, p, physics_inout%p)
	340	CALL pack_domain(info, pk, physics_inout%pk)
	341	CALL pack_domain(info, Temp, physics_inout%Temp)
	342	CALL pack_domain(info, ulon, physics_inout%ulon)
	343	CALL pack_domain(info, ulat, physics_inout%ulat)
	344	CALL pack_domain(info, q, physics_inout%q)
	345	ENDIF
	346	!$OMP BARRIER
[213]	347	END SUBROUTINE pack_physics
	348
[295]	349	SUBROUTINE unpack_physics(info, ps,temp, q, dulon, dulat)
[213]	350	USE theta2theta_rhodz_mod
	351	TYPE(t_pack_info) :: info
	352	REAL(rstd) :: ps(iim*jjm)
[295]	353	REAL(rstd) :: temp(iim*jjm,llm)
[213]	354	REAL(rstd) :: q(iim*jjm,llm,nqtot)
	355	REAL(rstd) :: dulon(iim*jjm,llm)
	356	REAL(rstd) :: dulat(iim*jjm,llm)
	357
	358	REAL(rstd) :: dq(iim*jjm,llm,nqtot)
	359	REAL(rstd) :: dTemp(iim*jjm,llm)
[481]	360
	361	!$OMP BARRIER
	362	IF (is_omp_level_master) THEN
	363	CALL unpack_domain(info, dulon, physics_inout%dulon)
	364	CALL unpack_domain(info, dulat, physics_inout%dulat)
	365	CALL unpack_domain(info, dq, physics_inout%dq)
	366	CALL unpack_domain(info, Temp, physics_inout%Temp)
	367	CALL unpack_domain(info, dTemp, physics_inout%dTemp)
	368	q = q + physics_inout%dt_phys * dq
	369	Temp = Temp + physics_inout%dt_phys * dTemp
	370	ENDIF
	371	!$OMP BARRIER
	372
[295]	373	! CALL compute_temperature2theta_rhodz(ps,Temp,theta_rhodz,0)
[213]	374	END SUBROUTINE unpack_physics
	375
	376	SUBROUTINE compute_update_velocity(dulon, dulat, ue)
	377	USE wind_mod
	378	REAL(rstd) :: dulon(iim*jjm,llm)
	379	REAL(rstd) :: dulat(iim*jjm,llm)
	380	REAL(rstd) :: ue(3iimjjm,llm)
[584]	381	REAL(rstd) :: duc(iim*jjm,llm,3)
[213]	382	REAL(rstd) :: dt2, due
	383	INTEGER :: i,j,ij,l
	384	! Reconstruct wind tendencies at edges and add to normal wind
	385	CALL compute_wind_centered_from_lonlat_compound(dulon,dulat,duc)
	386	dt2=.5*physics_inout%dt_phys
[481]	387	DO l=ll_begin,ll_end
[213]	388	DO j=jj_begin,jj_end
	389	DO i=ii_begin,ii_end
	390	ij=(j-1)*iim+i
[584]	391	due = sum( (duc(ij,l,:) + duc(ij+t_right,l,:))*ep_e(ij+u_right,:) )
[213]	392	ue(ij+u_right,l) = ue(ij+u_right,l) + dt2*due
	393
[584]	394	due = sum( (duc(ij,l,:) + duc(ij+t_lup,l,:))*ep_e(ij+u_lup,:) )
[213]	395	ue(ij+u_lup,l)=ue(ij+u_lup,l) + dt2*due
	396
[584]	397	due = sum( (duc(ij,l,:) + duc(ij+t_ldown,l,:))*ep_e(ij+u_ldown,:) )
[213]	398	ue(ij+u_ldown,l)=ue(ij+u_ldown,l) + dt2*due
	399	ENDDO
	400	ENDDO
	401	ENDDO
	402	END SUBROUTINE compute_update_velocity
	403
[81]	404	END MODULE physics_mod

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

Context Navigation

source: codes/icosagcm/devel/src/physics/physics.f90

Download in other formats: