Context Navigation

timestep_unstructured.F90 @ 833

Last change on this file since 833 was 833, checked in by dubos, 5 years ago
devel : fix pure-MPI build
File size: 12.0 KB

Rev	Line
[638]	1	MODULE timestep_unstructured_mod
	2	USE ISO_C_BINDING
	3	#ifdef CPP_USING_XIOS
	4	USE xios
	5	#endif
[833]	6	USE data_unstructured_mod
[638]	7	USE caldyn_unstructured_mod
[681]	8	USE transfer_unstructured_mod
[638]	9	IMPLICIT NONE
	10	PRIVATE
	11	SAVE
	12
	13	#ifdef CPP_USING_XIOS
	14	TYPE(xios_context) :: ctx_hdl
	15	#endif
	16
	17	CONTAINS
	18
[688]	19	#include "unstructured.h90"
[638]	20
	21	#define HASNAN(field) (ANY(.NOT.ABS(field)<1e20))
	22
[784]	23	SUBROUTINE scalar_laplacian(s,laps) BINDC(scalar_laplacian)
	24	FIELD_MASS :: s,laps ! IN, OUT
	25	FIELD_U :: grads ! BUF
	26	!$OMP PARALLEL NUM_THREADS(nb_threads)
	27	CALL update_halo(transfer_primal, s)
	28	CALL compute_scalar_laplacian(s,grads,laps)
	29	!$OMP END PARALLEL
	30	END SUBROUTINE scalar_laplacian
	31
[792]	32	SUBROUTINE curl_laplacian(u,lapu) BINDC(curl_laplacian)
	33	FIELD_U :: u,lapu ! IN, OUT
	34	FIELD_Z :: qv ! BUF
	35	!$OMP PARALLEL NUM_THREADS(nb_threads)
	36	CALL update_halo(transfer_edge, u)
	37	CALL compute_curl_laplacian(u,qv,lapu)
	38	!$OMP END PARALLEL
	39	END SUBROUTINE curl_laplacian
	40
[784]	41	!
[642]	42	SUBROUTINE caldyn_unstructured(tau, mass_col,rhodz,theta_rhodz,u,geopot,w, & ! IN : flow state
	43	theta,ps,pk,hflux,qv, & ! OUT : diags (except surface geopot : IN)
	44	dmass_col,drhodz,dtheta_rhodz,du_fast,du_slow, &
	45	dPhi_fast, dPhi_slow, dW_fast, dW_slow) BINDC(caldyn_unstructured) ! OUT : tendencies
[688]	46	NUM, VALUE :: tau
[645]	47	FIELD_MASS :: rhodz, drhodz, pk, berni ! IN, OUT, DIAG
	48	FIELD_THETA :: theta_rhodz, dtheta_rhodz, theta ! IN, OUT, DIAG
	49	FIELD_GEOPOT :: wflux, w, geopot, & ! DIAG, INOUT
	50	dPhi_fast, dPhi_slow, dW_fast, dW_slow ! OUT
[665]	51	FIELD_U :: u,du_fast,du_slow,hflux,qu ! INOUT,OUT,OUT,DIAG
[645]	52	FIELD_Z :: qv ! DIAG
	53	FIELD_PS :: ps,dmass_col,mass_col ! OUT,OUT,IN (if eta_mass) or OUT,UNUSED,UNUSED (if eta_lag)
[665]	54	! buffers for fields that need to be shared among OpenMP threads
	55	! vertical size is allocated as llm+1 even if only llm is needed
	56	FIELD_GEOPOT :: bufm1, bufm2, bufm3
	57	FIELD_UL :: bufu1, bufu2, bufu3,bufu4
	58
[688]	59	TIME :: time1,time2
[642]	60	INTEGER :: ij
	61
	62	time1=OMP_GET_WTIME()
	63
	64	IF(hydrostatic) THEN
	65
	66	!$OMP PARALLEL NUM_THREADS(nb_threads)
	67	!$OMP DO SCHEDULE(STATIC)
	68	DO ij=1,edge_num
	69	du_fast(:,ij)=0.
	70	du_slow(:,ij)=0.
	71	END DO
	72	!$OMP END DO
	73	CALL compute_theta(mass_col,rhodz,theta_rhodz, theta)
	74	CALL compute_geopot(rhodz,theta, ps,pk,geopot)
	75
	76	CALL compute_caldyn_fast(tau, pk,berni,theta,geopot, du_fast,u)
	77
	78	CALL compute_pvort_only(rhodz,u,qv,qu)
	79	CALL compute_caldyn_slow_hydro(rhodz,theta,u, berni,hflux,du_slow)
[665]	80	CALL compute_coriolis(hflux,theta,qu, bufu1, drhodz,dtheta_rhodz,du_slow)
[642]	81	IF(caldyn_eta == eta_mass) THEN
[665]	82	CALL caldyn_vert(drhodz,rhodz,theta,u, dmass_col,wflux,dtheta_rhodz,du_slow,bufu1)
[642]	83	END IF
	84	!$OMP END PARALLEL
	85
	86	ELSE ! NH
	87	DO ij=1,edge_num
	88	du_fast(:,ij)=0.
	89	du_slow(:,ij)=0.
	90	END DO
	91	DO ij=1,primal_num
	92	wflux(1,ij)=0.
	93	wflux(llm+1,ij)=0.
	94	END DO
	95	CALL compute_theta(mass_col,rhodz,theta_rhodz, theta)
[665]	96	CALL compute_caldyn_solver(tau,rhodz,theta, bufm1,bufm2,bufm3, pk,geopot,W,dPhi_fast,dW_fast,du_fast)
[642]	97	CALL compute_caldyn_fast(tau, pk,berni,theta,geopot, du_fast,u)
	98	CALL compute_pvort_only(rhodz,u,qv,qu)
[665]	99	CALL compute_caldyn_slow_NH(u,rhodz,geopot,W, bufm1,bufm2,bufm3,bufu1,bufu2,bufu3,bufu4, hflux,du_slow,dPhi_slow,dW_slow)
	100	CALL compute_coriolis(hflux,theta,qu, bufu1, drhodz,dtheta_rhodz,du_slow)
[642]	101	IF(caldyn_eta == eta_mass) THEN
[665]	102	CALL caldyn_vert(drhodz,rhodz,theta,u, dmass_col,wflux,dtheta_rhodz,du_slow,bufu1)
	103	CALL compute_caldyn_vert_NH(rhodz,geopot,W,wflux, bufm1,bufm2,bufm3, du_slow,dPhi_slow,dW_slow)
[642]	104	END IF
[638]	105	END IF
[642]	106
	107	time2=OMP_GET_WTIME()
	108	IF(time2>time1) elapsed = elapsed + time2-time1
[652]	109
	110	CALL print_trace()
[642]	111	END SUBROUTINE caldyn_unstructured
	112	!
	113	!----------------------------- Time stepping -------------------------------
	114
	115	!
	116	SUBROUTINE ARK_step(nstep, mass_col,rhodz,theta_rhodz,u,geopot,w, & ! INOUT : flow state
[638]	117	theta,ps,pk,hflux,qv, & ! OUT : diags (except surface geopot : IN)
	118	dmass_col,drhodz,dtheta_rhodz,du_fast,du_slow, &
	119	dPhi_fast, dPhi_slow, dW_fast, dW_slow) BINDC(ARK_step) ! OUT : tendencies
[642]	120	INTEGER, VALUE :: nstep ! advance by nstep time steps
[638]	121	FIELD_PS :: mass_col, ps ! OUT,IN (if eta_mass) or OUT,UNUSED (if eta_lag)
	122	FIELD_MASS :: rhodz, pk, berni ! IN, DIAG
	123	FIELD_THETA :: theta_rhodz, theta ! IN, DIAG
[665]	124	FIELD_U :: u,hflux,qu ! INOUT,DIAG*2
	125	FIELD_GEOPOT :: geopot, w, wflux ! IN, INOUT, DIAG
[638]	126	FIELD_Z :: qv ! DIAG
[688]	127	NUM2( primal_num, max_nb_stage) :: dmass_col ! OUT
	128	NUM3(llm, primal_num, max_nb_stage) :: drhodz ! OUT
	129	NUM4(llm, primal_num, nqdyn, max_nb_stage) :: dtheta_rhodz ! OUT
	130	NUM3(llm, edge_num, max_nb_stage) :: du_fast,du_slow ! OUT
	131	NUM3(llm+1, primal_num, max_nb_stage) :: &
[638]	132	dPhi_fast, dPhi_slow, dW_fast, dW_slow ! OUT
[665]	133	! buffers for fields that need to be shared among OpenMP threads
	134	! vertical size is allocated as llm+1 even if only llm is needed
	135	FIELD_GEOPOT :: bufm1, bufm2, bufm3
	136	FIELD_UL :: bufu1, bufu2, bufu3,bufu4
[688]	137	TIME :: time1,time2
[683]	138	INTEGER :: step, stage, iq, ij
[642]	139
[638]	140	time1=OMP_GET_WTIME()
[642]	141
[638]	142	!$OMP PARALLEL NUM_THREADS(nb_threads)
	143
[677]	144	!$OMP DO SCHEDULE(STATIC)
	145	DO ij=1,primal_num
	146	wflux(1,ij)=0.
	147	wflux(llm+1,ij)=0.
	148	END DO
	149	!$OMP END DO
	150
[642]	151	DO step=1, nstep
	152	DO stage=1, nb_stage
[638]	153
[681]	154	CALL update_halo(transfer_primal, rhodz)
[683]	155	DO iq=1, nqdyn
	156	CALL update_halo(transfer_primal, theta_rhodz(:,:,iq))
	157	END DO
[681]	158
[642]	159	IF(hydrostatic) THEN
	160
	161	!$OMP DO SCHEDULE(STATIC)
	162	DO ij=1,edge_num
	163	du_fast(:,ij,stage)=0.
	164	du_slow(:,ij,stage)=0.
	165	END DO
[775]	166	!$OMP END DO
[642]	167
	168	CALL compute_theta(mass_col,rhodz,theta_rhodz, theta)
	169	CALL compute_geopot(rhodz,theta, ps,pk,geopot)
	170
	171	CALL compute_caldyn_fast(tauj(stage), pk,berni,theta,geopot, du_fast(:,:,stage), u)
[681]	172	CALL update_halo(transfer_edge, u)
[642]	173	CALL compute_pvort_only(rhodz,u,qv,qu)
[681]	174	CALL update_halo(transfer_edge, qu)
[775]	175
[642]	176	CALL compute_caldyn_slow_hydro(rhodz,theta,u, berni,hflux,du_slow(:,:,stage))
[683]	177	CALL compute_coriolis(hflux,theta,qu, bufu1, drhodz(:,:,stage), &
	178	dtheta_rhodz(:,:,:,stage),du_slow(:,:,stage))
[642]	179	IF(caldyn_eta == eta_mass) THEN
	180	CALL caldyn_vert(drhodz(:,:,stage),rhodz,theta,u, &
[665]	181	dmass_col(:,stage),wflux,dtheta_rhodz(:,:,:,stage),du_slow(:,:,stage),bufu1)
[642]	182	END IF
	183
	184	ELSE ! NH
[683]	185	CALL update_halo(transfer_primal, geopot)
	186	CALL update_halo(transfer_primal, W)
	187
[642]	188	CALL compute_theta(mass_col,rhodz,theta_rhodz, theta)
[665]	189	CALL compute_caldyn_solver(tauj(stage),rhodz,theta, bufm1,bufm2,bufm3, pk,geopot,W, &
[642]	190	dPhi_fast(:,:,stage), dW_fast(:,:,stage), du_fast(:,:,stage))
[683]	191
[642]	192	CALL compute_caldyn_fast(tauj(stage), pk,berni,theta,geopot, du_fast(:,:,stage),u)
[683]	193	CALL update_halo(transfer_edge, u)
[642]	194	CALL compute_pvort_only(rhodz,u,qv,qu)
[683]	195	CALL update_halo(transfer_edge, qu)
	196
	197	CALL compute_caldyn_slow_NH(u,rhodz,geopot,W, bufm1,bufm2,bufm3, &
	198	bufu1,bufu2,bufu3,bufu4, hflux, &
[665]	199	du_slow(:,:,stage), dPhi_slow(:,:,stage), dW_slow(:,:,stage))
[683]	200	CALL compute_coriolis(hflux,theta,qu, bufu1, &
	201	drhodz(:,:,stage),dtheta_rhodz(:,:,:,stage), du_slow(:,:,stage))
[642]	202	IF(caldyn_eta == eta_mass) THEN
	203	CALL caldyn_vert(drhodz(:,:,stage),rhodz,theta,u, &
[665]	204	dmass_col(:,stage),wflux,dtheta_rhodz(:,:,:,stage), du_slow(:,:,stage), bufu1)
	205	CALL compute_caldyn_vert_NH(rhodz,geopot,W,wflux, bufm1,bufm2,bufm3, &
	206	du_slow(:,:,stage), dPhi_slow(:,:,stage), dW_slow(:,:,stage) )
[642]	207	END IF
	208	END IF ! NH
[638]	209
[642]	210	! FIXME : mass_col is computed from rhodz
	211	! so that the DOFs are the same whatever caldyn_eta
	212	! in DYNAMICO mass_col is prognosed rather than rhodz
	213
	214	#define UPDATE(clj,field,dfield) update(stage,SIZE(field),clj,field,dfield)
	215
	216	CALL UPDATE(cslj, rhodz, drhodz)
	217	CALL UPDATE(cslj, theta_rhodz, dtheta_rhodz)
	218	CALL UPDATE(cslj, u, du_slow)
	219	CALL UPDATE(cflj, u, du_fast)
	220
	221	IF(.NOT.hydrostatic) THEN
	222	CALL UPDATE(cslj, geopot, dPhi_slow)
	223	CALL UPDATE(cflj, geopot, dPhi_fast)
	224	CALL UPDATE(cslj, W, dW_slow)
	225	CALL UPDATE(cflj, W, dW_fast)
[638]	226	END IF
	227
[642]	228	END DO
[638]	229	END DO
[642]	230	!$OMP END PARALLEL
[638]	231
[642]	232	time2=OMP_GET_WTIME()
	233	IF(time2>time1) elapsed = elapsed + time2-time1
[652]	234
	235	CALL print_trace()
[642]	236
[638]	237	END SUBROUTINE ARK_step
[642]	238
	239
	240	SUBROUTINE update(j,sz,clj,field,dfield)
	241	INTEGER :: j, sz ! stage in ARK scheme, field size
[688]	242	NUM2(max_nb_stage,max_nb_stage) :: clj ! modified Butcher tableau
	243	NUM1(sz) :: field
	244	NUM2(sz, max_nb_stage) :: dfield
[642]	245	!
	246	INTEGER :: l, ij
[658]	247	CALL enter_trace(id_update, 8sz(j+1) )
[642]	248	! PRINT *, clj(1:j,j)
	249	SELECT CASE(j)
	250	CASE(1)
	251	!$OMP DO SCHEDULE(static)
	252	DO ij=1,sz
	253	field(ij) = field(ij) &
	254	+ clj(1,j)*dfield(ij,1)
	255	END DO
	256	CASE(2)
	257	!$OMP DO SCHEDULE(static)
	258	DO ij=1,sz
	259	field(ij) = field(ij) &
	260	+ clj(1,j)*dfield(ij,1) &
	261	+ clj(2,j)*dfield(ij,2)
	262	END DO
	263	CASE(3)
	264	!$OMP DO SCHEDULE(static)
	265	DO ij=1,sz
	266	field(ij) = field(ij) &
	267	+ clj(1,j)*dfield(ij,1) &
	268	+ clj(2,j)*dfield(ij,2) &
	269	+ clj(3,j)*dfield(ij,3)
	270
	271	END DO
	272	CASE(4)
	273	!$OMP DO SCHEDULE(static)
	274	DO ij=1,sz
	275	field(ij) = field(ij) &
	276	+ clj(1,j)*dfield(ij,1) &
	277	+ clj(2,j)*dfield(ij,2) &
	278	+ clj(3,j)*dfield(ij,3) &
	279	+ clj(4,j)*dfield(ij,4)
	280	END DO
	281	CASE(5)
	282	!$OMP DO SCHEDULE(static)
	283	DO ij=1,sz
	284	field(ij) = field(ij) &
	285	+ clj(1,j)*dfield(ij,1) &
	286	+ clj(2,j)*dfield(ij,2) &
	287	+ clj(3,j)*dfield(ij,3) &
	288	+ clj(4,j)*dfield(ij,4) &
	289	+ clj(5,j)*dfield(ij,5)
	290	END DO
	291	END SELECT
[658]	292	CALL exit_trace()
[642]	293	END SUBROUTINE update
	294
	295	!---------------------------------------------- XIOS ----------------------------------------
[638]	296
	297	#ifdef CPP_USING_XIOS
[642]	298
[638]	299	SUBROUTINE setup_xios() BINDC(setup_xios)
	300	! MPI_INIT / MPI_finalize are assumed to be called BEFORE/AFTER this routine
	301	! This is the case when calling from Python after importing mpi4py
	302	INTEGER :: ierr, mpi_threading_mode
	303	PRINT *, 'Initialize XIOS and obtain our communicator'
	304	CALL xios_initialize("icosagcm",return_comm=comm_icosa)
	305	CALL xios_context_initialize("icosagcm",comm_icosa)
	306	CALL xios_get_handle("icosagcm",ctx_hdl)
	307	CALL xios_set_current_context(ctx_hdl)
	308	END SUBROUTINE setup_xios
[683]	309
[638]	310	SUBROUTINE call_xios_set_timestep(dt) BINDC(xios_set_timestep)
	311	DBL, VALUE :: dt
	312	TYPE(xios_duration) :: dtime
	313	dtime%second=dt
	314	CALL xios_set_timestep(dtime)
	315	END SUBROUTINE call_xios_set_timestep
	316
	317	SUBROUTINE call_xios_update_calendar(step) BINDC(xios_update_calendar)
	318	INDEX, value :: step
	319	CALL xios_update_calendar(step)
	320	END SUBROUTINE call_xios_update_calendar
	321
	322	#endif
	323
	324	END MODULE timestep_unstructured_mod

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

Context Navigation

source: codes/icosagcm/devel/src/unstructured/timestep_unstructured.F90 @ 833

Download in other formats: