Context Navigation

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

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

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

Download in other formats: