Context Navigation

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

Last change on this file since 775 was 775, checked in by dubos, 6 years ago
devel/unstructured : fix missing OMP END DO
File size: 11.5 KB

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

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

Download in other formats: