Context Navigation

source: codes/icosagcm/trunk/src/timeloop_gcm.f90 @ 167

Last change on this file since 167 was 167, checked in by dubos, 11 years ago
Multi-layer Saint-Venant / Ripa equations - tested with Williamson91.5 & JW06
File size: 19.8 KB

Rev	Line
[12]	1	MODULE timeloop_gcm_mod
[151]	2	USE transfert_mod
	3	USE icosa
[133]	4	PRIVATE
[12]	5
[151]	6	PUBLIC :: init_timeloop, timeloop
[133]	7
	8	INTEGER, PARAMETER :: euler=1, rk4=2, mlf=3
[148]	9	INTEGER :: itau_sync=10
[133]	10
[162]	11	TYPE(t_message) :: req_ps0, req_mass0, req_theta_rhodz0, req_u0, req_q0
[151]	12
	13	TYPE(t_field),POINTER :: f_q(:)
[159]	14	TYPE(t_field),POINTER :: f_rhodz(:), f_mass(:), f_massm1(:), f_massm2(:), f_dmass(:)
	15	TYPE(t_field),POINTER :: f_phis(:), f_ps(:),f_psm1(:), f_psm2(:), f_dps(:)
	16	TYPE(t_field),POINTER :: f_u(:),f_um1(:),f_um2(:), f_du(:)
	17	TYPE(t_field),POINTER :: f_theta_rhodz(:),f_theta_rhodzm1(:),f_theta_rhodzm2(:), f_dtheta_rhodz(:)
[151]	18	TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:), f_hfluxt(:), f_wfluxt(:)
	19
	20	INTEGER :: nb_stage, matsuno_period, scheme
	21
	22	REAL(rstd),SAVE :: jD_cur, jH_cur
	23	REAL(rstd),SAVE :: start_time
	24
[12]	25	CONTAINS
	26
[151]	27	SUBROUTINE init_timeloop
[19]	28	USE icosa
[15]	29	USE dissip_gcm_mod
[17]	30	USE caldyn_mod
[12]	31	USE etat0_mod
[159]	32	USE disvert_mod
[17]	33	USE guided_mod
	34	USE advect_tracer_mod
[81]	35	USE physics_mod
[131]	36	USE mpipara
[151]	37	USE omp_para
[145]	38	USE trace
[148]	39	USE transfert_mod
[151]	40	USE check_conserve_mod
	41	USE ioipsl
[12]	42	IMPLICIT NONE
	43
[159]	44	CHARACTER(len=255) :: def
[17]	45
[149]	46	!----------------------------------------------------
	47	IF (TRIM(time_style)=='lmd') Then
	48
	49	day_step=180
	50	CALL getin('day_step',day_step)
	51
	52	ndays=1
	53	CALL getin('ndays',ndays)
	54
	55	dt = daysec/REAL(day_step)
	56	itaumax = ndays*day_step
	57
	58	calend = 'earth_360d'
	59	CALL getin('calend', calend)
	60
	61	day_ini = 0
	62	CALL getin('day_ini',day_ini)
	63
	64	day_end = 0
	65	CALL getin('day_end',day_end)
	66
	67	annee_ref = 1998
	68	CALL getin('annee_ref',annee_ref)
	69
	70	start_time = 0
	71	CALL getin('start_time',start_time)
	72
	73	write_period=0
	74	CALL getin('write_period',write_period)
	75
	76	write_period=write_period/scale_factor
	77	itau_out=FLOOR(write_period/dt)
	78
	79	PRINT *, 'Output frequency (scaled) set to ',write_period, ' : itau_out = ',itau_out
	80
	81	mois = 1 ; heure = 0.
	82	call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
	83	jH_ref = jD_ref - int(jD_ref)
	84	jD_ref = int(jD_ref)
	85
	86	CALL ioconf_startdate(INT(jD_ref),jH_ref)
	87	write(,)'annee_ref, mois, day_ref, heure, jD_ref'
	88	write(,)annee_ref, mois, day_ref, heure, jD_ref
	89	write(,)"ndays,day_step,itaumax,dt======>"
	90	write(,)ndays,day_step,itaumax,dt
	91	call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
	92	write(,)'jD_ref+jH_ref,an, mois, jour, heure'
	93	write(,)jD_ref+jH_ref,an, mois, jour, heure
	94	day_end = day_ini + ndays
	95	END IF
	96	!----------------------------------------------------
	97
[129]	98
[159]	99	! Time-independant orography
	100	CALL allocate_field(f_phis,field_t,type_real,name='phis')
	101	! Trends
	102	CALL allocate_field(f_du,field_u,type_real,llm,name='du')
	103	CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm,name='dtheta_rhodz')
	104	! Model state at current time step (RK/MLF/Euler)
	105	CALL allocate_field(f_ps,field_t,type_real, name='ps')
	106	CALL allocate_field(f_mass,field_t,type_real,llm,name='mass')
	107	CALL allocate_field(f_u,field_u,type_real,llm,name='u')
	108	CALL allocate_field(f_theta_rhodz,field_t,type_real,llm,name='theta_rhodz')
	109	! Model state at previous time step (RK/MLF)
	110	CALL allocate_field(f_um1,field_u,type_real,llm,name='um1')
	111	CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm,name='theta_rhodzm1')
	112	! Tracers
	113	CALL allocate_field(f_q,field_t,type_real,llm,nqtot)
	114	CALL allocate_field(f_rhodz,field_t,type_real,llm,name='rhodz')
	115	! Mass fluxes
	116	CALL allocate_field(f_hflux,field_u,type_real,llm) ! instantaneous mass fluxes
	117	CALL allocate_field(f_hfluxt,field_u,type_real,llm) ! mass "fluxes" accumulated in time
	118	CALL allocate_field(f_wflux,field_t,type_real,llm+1) ! vertical mass fluxes
[162]	119	CALL allocate_field(f_dmass,field_t,type_real,llm, name='dmass')
[151]	120
[159]	121	IF(caldyn_eta == eta_mass) THEN ! eta = mass coordinate (default)
	122	CALL allocate_field(f_dps,field_t,type_real,name='dps')
	123	CALL allocate_field(f_psm1,field_t,type_real,name='psm1')
	124	CALL allocate_field(f_wfluxt,field_t,type_real,llm+1,name='wfluxt')
	125	! the following are unused but must point to something
[162]	126	! f_massm1 => f_mass
[159]	127	ELSE ! eta = Lagrangian vertical coordinate
[162]	128	CALL allocate_field(f_massm1,field_t,type_real,llm, name='massm1')
[159]	129	! the following are unused but must point to something
	130	f_wfluxt => f_wflux
	131	f_dps => f_phis
	132	f_psm1 => f_phis
	133	END IF
[151]	134
[159]	135	def='runge_kutta'
	136	CALL getin('scheme',def)
	137
	138	SELECT CASE (TRIM(def))
[151]	139	CASE('euler')
	140	scheme=euler
	141	nb_stage=1
	142	CASE ('runge_kutta')
	143	scheme=rk4
	144	nb_stage=4
	145	CASE ('leapfrog_matsuno')
	146	scheme=mlf
	147	matsuno_period=5
	148	CALL getin('matsuno_period',matsuno_period)
	149	nb_stage=matsuno_period+1
[129]	150	! Model state 2 time steps ago (MLF)
[151]	151	CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm)
	152	CALL allocate_field(f_um2,field_u,type_real,llm)
[159]	153	IF(caldyn_eta == eta_mass) THEN ! eta = mass coordinate (default)
	154	CALL allocate_field(f_psm2,field_t,type_real)
	155	! the following are unused but must point to something
	156	f_massm2 => f_mass
	157	ELSE ! eta = Lagrangian vertical coordinate
	158	CALL allocate_field(f_massm2,field_t,type_real,llm)
	159	! the following are unused but must point to something
	160	f_psm2 => f_phis
	161	END IF
	162
[151]	163	CASE default
[159]	164	PRINT*,'Bad selector for variable scheme : <', TRIM(def), &
[151]	165	' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>'
	166	STOP
	167	END SELECT
	168
	169
	170	CALL init_dissip
	171	CALL init_caldyn
	172	CALL init_guided
	173	CALL init_advect_tracer
	174	CALL init_check_conserve
	175	CALL init_physics
[159]	176	CALL etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_q)
[151]	177
	178	CALL transfert_request(f_phis,req_i0)
	179	CALL transfert_request(f_phis,req_i1)
	180	CALL writefield("phis",f_phis,once=.TRUE.)
	181
	182	CALL init_message(f_ps,req_i0,req_ps0)
[162]	183	CALL init_message(f_mass,req_i0,req_mass0)
[151]	184	CALL init_message(f_theta_rhodz,req_i0,req_theta_rhodz0)
	185	CALL init_message(f_u,req_e0_vect,req_u0)
	186	CALL init_message(f_q,req_i0,req_q0)
	187
	188	END SUBROUTINE init_timeloop
[12]	189
[151]	190	SUBROUTINE timeloop
	191	USE icosa
	192	USE dissip_gcm_mod
[159]	193	USE disvert_mod
[151]	194	USE caldyn_mod
[162]	195	USE caldyn_gcm_mod, ONLY : req_ps, req_mass
[151]	196	USE etat0_mod
	197	USE guided_mod
	198	USE advect_tracer_mod
	199	USE physics_mod
	200	USE mpipara
	201	USE omp_para
	202	USE trace
	203	USE transfert_mod
	204	USE check_conserve_mod
	205	IMPLICIT NONE
	206	REAL(rstd),POINTER :: q(:,:,:)
[159]	207	REAL(rstd),POINTER :: phis(:), ps(:) ,psm1(:), psm2(:), dps(:)
	208	REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:), du(:,:)
	209	REAL(rstd),POINTER :: rhodz(:,:), mass(:,:), massm1(:,:), massm2(:,:), dmass(:,:)
	210	REAL(rstd),POINTER :: theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:), dtheta_rhodz(:,:)
[151]	211	REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:)
[12]	212
[151]	213	INTEGER :: ind
	214	INTEGER :: it,i,j,n, stage
	215	LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time
	216	LOGICAL, PARAMETER :: check=.FALSE.
[50]	217
[159]	218	CALL caldyn_BC(f_phis, f_wflux) ! set constant values in first/last interfaces
[157]	219
[151]	220	!$OMP BARRIER
[133]	221	DO ind=1,ndomain
	222	CALL swap_dimensions(ind)
	223	CALL swap_geometry(ind)
[162]	224	rhodz=f_rhodz(ind); mass=f_mass(ind); ps=f_ps(ind)
	225	IF(caldyn_eta==eta_mass) THEN
	226	CALL compute_rhodz(.TRUE., ps, rhodz) ! save rhodz for transport scheme before dynamics update ps
	227	ELSE
	228	rhodz(:,:)=mass(:,:)
	229	END IF
[133]	230	END DO
[138]	231	fluxt_zero=.TRUE.
[132]	232
[12]	233	DO it=0,itaumax
[148]	234	IF (MOD(it,itau_sync)==0) THEN
[151]	235	CALL send_message(f_ps,req_ps0)
[162]	236	CALL send_message(f_mass,req_mass0)
[151]	237	CALL send_message(f_theta_rhodz,req_theta_rhodz0)
	238	CALL send_message(f_u,req_u0)
	239	CALL send_message(f_q,req_q0)
	240	CALL wait_message(req_ps0)
[162]	241	CALL wait_message(req_mass0)
[151]	242	CALL wait_message(req_theta_rhodz0)
	243	CALL wait_message(req_u0)
	244	CALL wait_message(req_q0)
[148]	245	ENDIF
	246
[151]	247	! IF (is_mpi_root) PRINT ,"It No :",It," t :",dtIt
[63]	248	IF (mod(it,itau_out)==0 ) THEN
[151]	249	CALL writefield("q",f_q)
[81]	250	CALL update_time_counter(dt*it)
[151]	251	CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)
[63]	252	ENDIF
[151]	253
	254	CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
[129]	255
	256	DO stage=1,nb_stage
	257	CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), &
[159]	258	f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q, &
[162]	259	f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du)
[133]	260	SELECT CASE (scheme)
	261	CASE(euler)
	262	CALL euler_scheme(.TRUE.)
	263	CASE (rk4)
[129]	264	CALL rk_scheme(stage)
[133]	265	CASE (mlf)
[129]	266	CALL leapfrog_matsuno_scheme(stage)
[133]	267	CASE DEFAULT
[129]	268	STOP
	269	END SELECT
	270	END DO
[130]	271
[148]	272	IF (MOD(it+1,itau_dissip)==0) THEN
[167]	273	IF(caldyn_eta==eta_mass) THEN
	274	DO ind=1,ndomain
	275	CALL swap_dimensions(ind)
	276	CALL swap_geometry(ind)
	277	mass=f_mass(ind); ps=f_ps(ind);
	278	CALL compute_rhodz(.TRUE., ps, mass)
	279	END DO
	280	ENDIF
	281	CALL dissip(f_u,f_du,f_mass,f_phis, f_theta_rhodz,f_dtheta_rhodz)
	282	CALL euler_scheme(.FALSE.) ! update only u, theta
	283	END IF
[130]	284
[133]	285	IF(MOD(it+1,itau_adv)==0) THEN
[138]	286
[135]	287	CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz) ! update q and rhodz after RK step
[134]	288	fluxt_zero=.TRUE.
[138]	289
	290	! FIXME : check that rhodz is consistent with ps
[148]	291	IF (check) THEN
	292	DO ind=1,ndomain
	293	CALL swap_dimensions(ind)
	294	CALL swap_geometry(ind)
[151]	295	rhodz=f_rhodz(ind); ps=f_ps(ind);
[148]	296	CALL compute_rhodz(.FALSE., ps, rhodz)
	297	END DO
	298	ENDIF
[151]	299
[133]	300	END IF
[151]	301
	302
	303
[149]	304	!----------------------------------------------------
[151]	305	! jD_cur = jD_ref + day_ini - day_ref + it/day_step
	306	! jH_cur = jH_ref + start_time + mod(it,day_step)/float(day_step)
	307	! jD_cur = jD_cur + int(jH_cur)
	308	! jH_cur = jH_cur - int(jH_cur)
	309	! CALL physics(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
	310
[124]	311	! CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
[129]	312	ENDDO
[151]	313
[129]	314
[12]	315	CONTAINS
	316
[130]	317	SUBROUTINE Euler_scheme(with_dps)
[12]	318	IMPLICIT NONE
[130]	319	LOGICAL :: with_dps
	320	INTEGER :: ind
[148]	321	INTEGER :: i,j,ij,l
[145]	322	CALL trace_start("Euler_scheme")
	323
[130]	324	DO ind=1,ndomain
[138]	325	CALL swap_dimensions(ind)
	326	CALL swap_geometry(ind)
[148]	327
[162]	328	IF(with_dps) THEN ! update ps/mass
	329	IF(caldyn_eta==eta_mass) THEN ! update ps
	330	ps=f_ps(ind) ; dps=f_dps(ind) ;
	331	IF (omp_first) THEN
	332	DO j=jj_begin,jj_end
	333	DO i=ii_begin,ii_end
	334	ij=(j-1)*iim+i
	335	ps(ij)=ps(ij)+dt*dps(ij)
	336	ENDDO
	337	ENDDO
	338	ENDIF
	339	ELSE ! update mass
	340	mass=f_mass(ind) ; dmass=f_dmass(ind) ;
	341	DO l=1,llm
	342	DO j=jj_begin,jj_end
	343	DO i=ii_begin,ii_end
	344	ij=(j-1)*iim+i
	345	mass(ij,l)=mass(ij,l)+dt*dmass(ij,l)
	346	ENDDO
	347	ENDDO
	348	END DO
	349	END IF
	350
	351	hflux=f_hflux(ind); hfluxt=f_hfluxt(ind)
	352	wflux=f_wflux(ind); wfluxt=f_wfluxt(ind)
	353	CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero(ind))
	354	END IF ! update ps/mass
[148]	355
[130]	356	u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
	357	du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
[148]	358
[151]	359	DO l=ll_begin,ll_end
[148]	360	DO j=jj_begin,jj_end
	361	DO i=ii_begin,ii_end
	362	ij=(j-1)*iim+i
	363	u(ij+u_right,l)=u(ij+u_right,l)+dt*du(ij+u_right,l)
	364	u(ij+u_lup,l)=u(ij+u_lup,l)+dt*du(ij+u_lup,l)
	365	u(ij+u_ldown,l)=u(ij+u_ldown,l)+dt*du(ij+u_ldown,l)
	366	theta_rhodz(ij,l)=theta_rhodz(ij,l)+dt*dtheta_rhodz(ij,l)
	367	ENDDO
	368	ENDDO
	369	ENDDO
[130]	370	ENDDO
[133]	371
[145]	372	CALL trace_end("Euler_scheme")
	373
[12]	374	END SUBROUTINE Euler_scheme
[120]	375
[129]	376	SUBROUTINE RK_scheme(stage)
[120]	377	IMPLICIT NONE
	378	INTEGER :: ind, stage
[129]	379	REAL(rstd), DIMENSION(4), PARAMETER :: coef = (/ .25, 1./3., .5, 1. /)
[120]	380	REAL(rstd) :: tau
[148]	381	INTEGER :: i,j,ij,l
[145]	382
	383	CALL trace_start("RK_scheme")
[120]	384
	385	tau = dt*coef(stage)
[151]	386
[159]	387	! if mass coordinate, deal with ps first on one core
	388	IF(caldyn_eta==eta_mass) THEN
	389	IF (omp_first) THEN
[162]	390
[159]	391	DO ind=1,ndomain
	392	CALL swap_dimensions(ind)
	393	CALL swap_geometry(ind)
[162]	394	ps=f_ps(ind) ; psm1=f_psm1(ind) ; dps=f_dps(ind)
[159]	395
	396	IF (stage==1) THEN ! first stage : save model state in XXm1
	397	DO j=jj_begin,jj_end
	398	DO i=ii_begin,ii_end
	399	ij=(j-1)*iim+i
	400	psm1(ij)=ps(ij)
	401	ENDDO
	402	ENDDO
	403	ENDIF
	404
	405	! updates are of the form x1 := x0 + tau*f(x1)
	406	DO j=jj_begin,jj_end
	407	DO i=ii_begin,ii_end
	408	ij=(j-1)*iim+i
	409	ps(ij)=psm1(ij)+tau*dps(ij)
	410	ENDDO
[151]	411	ENDDO
[159]	412	ENDDO
	413	ENDIF
	414	CALL send_message(f_ps,req_ps)
[162]	415
	416	ELSE ! Lagrangian coordinate, deal with mass
	417	DO ind=1,ndomain
	418	CALL swap_dimensions(ind)
	419	CALL swap_geometry(ind)
	420	mass=f_mass(ind); dmass=f_dmass(ind); massm1=f_massm1(ind)
	421
	422	IF (stage==1) THEN ! first stage : save model state in XXm1
	423	DO l=ll_begin,ll_end
	424	DO j=jj_begin,jj_end
	425	DO i=ii_begin,ii_end
	426	ij=(j-1)*iim+i
	427	massm1(ij,l)=mass(ij,l)
	428	ENDDO
	429	ENDDO
	430	ENDDO
	431	END IF
	432
	433	! updates are of the form x1 := x0 + tau*f(x1)
	434	DO l=ll_begin,ll_end
	435	DO j=jj_begin,jj_end
	436	DO i=ii_begin,ii_end
	437	ij=(j-1)*iim+i
	438	mass(ij,l)=massm1(ij,l)+tau*dmass(ij,l)
	439	ENDDO
	440	ENDDO
	441	ENDDO
	442	END DO
	443	CALL send_message(f_mass,req_mass)
	444
[159]	445	END IF
[151]	446
[159]	447	! now deal with other prognostic variables
[120]	448	DO ind=1,ndomain
[138]	449	CALL swap_dimensions(ind)
	450	CALL swap_geometry(ind)
[162]	451	u=f_u(ind) ; du=f_du(ind) ; um1=f_um1(ind)
	452	theta_rhodz=f_theta_rhodz(ind)
	453	theta_rhodzm1=f_theta_rhodzm1(ind)
	454	dtheta_rhodz=f_dtheta_rhodz(ind)
[129]	455
	456	IF (stage==1) THEN ! first stage : save model state in XXm1
[159]	457	DO l=ll_begin,ll_end
[148]	458	DO j=jj_begin,jj_end
	459	DO i=ii_begin,ii_end
	460	ij=(j-1)*iim+i
	461	um1(ij+u_right,l)=u(ij+u_right,l)
	462	um1(ij+u_lup,l)=u(ij+u_lup,l)
	463	um1(ij+u_ldown,l)=u(ij+u_ldown,l)
	464	theta_rhodzm1(ij,l)=theta_rhodz(ij,l)
	465	ENDDO
	466	ENDDO
	467	ENDDO
[162]	468	END IF
[148]	469
[151]	470	DO l=ll_begin,ll_end
[148]	471	DO j=jj_begin,jj_end
	472	DO i=ii_begin,ii_end
	473	ij=(j-1)*iim+i
	474	u(ij+u_right,l)=um1(ij+u_right,l)+tau*du(ij+u_right,l)
	475	u(ij+u_lup,l)=um1(ij+u_lup,l)+tau*du(ij+u_lup,l)
	476	u(ij+u_ldown,l)=um1(ij+u_ldown,l)+tau*du(ij+u_ldown,l)
	477	theta_rhodz(ij,l)=theta_rhodzm1(ij,l)+tau*dtheta_rhodz(ij,l)
	478	ENDDO
	479	ENDDO
	480	ENDDO
[162]	481
[133]	482	IF(stage==nb_stage) THEN ! accumulate mass fluxes at last stage
	483	hflux=f_hflux(ind); hfluxt=f_hfluxt(ind)
[138]	484	wflux=f_wflux(ind); wfluxt=f_wfluxt(ind)
	485	CALL accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, dt,fluxt_zero(ind))
[133]	486	END IF
[120]	487	END DO
[145]	488
	489	CALL trace_end("RK_scheme")
	490
[120]	491	END SUBROUTINE RK_scheme
	492
[129]	493	SUBROUTINE leapfrog_matsuno_scheme(stage)
[12]	494	IMPLICIT NONE
[129]	495	INTEGER :: ind, stage
	496	REAL :: tau
[145]	497
	498	CALL trace_start("leapfrog_matsuno_scheme")
	499
	500	tau = dt/nb_stage
[12]	501	DO ind=1,ndomain
[138]	502	CALL swap_dimensions(ind)
	503	CALL swap_geometry(ind)
	504
[12]	505	ps=f_ps(ind) ; u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
	506	psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
	507	psm2=f_psm2(ind) ; um2=f_um2(ind) ; theta_rhodzm2=f_theta_rhodzm2(ind)
	508	dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
	509
	510
[129]	511	IF (stage==1) THEN
[12]	512	psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:)
[129]	513	ps(:)=psm1(:)+tau*dps(:)
	514	u(:,:)=um1(:,:)+tau*du(:,:)
	515	theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
[12]	516
[129]	517	ELSE IF (stage==2) THEN
[12]	518
[129]	519	ps(:)=psm1(:)+tau*dps(:)
	520	u(:,:)=um1(:,:)+tau*du(:,:)
	521	theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
[12]	522
	523	psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:)
	524	psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:)
	525
	526	ELSE
	527
[129]	528	ps(:)=psm2(:)+2taudps(:)
	529	u(:,:)=um2(:,:)+2taudu(:,:)
	530	theta_rhodz(:,:)=theta_rhodzm2(:,:)+2taudtheta_rhodz(:,:)
[12]	531
	532	psm2(:)=psm1(:) ; theta_rhodzm2(:,:)=theta_rhodzm1(:,:) ; um2(:,:)=um1(:,:)
	533	psm1(:)=ps(:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:) ; um1(:,:)=u(:,:)
	534
	535	ENDIF
	536
	537	ENDDO
[145]	538	CALL trace_end("leapfrog_matsuno_scheme")
[12]	539
	540	END SUBROUTINE leapfrog_matsuno_scheme
	541
[50]	542	END SUBROUTINE timeloop
[133]	543
[159]	544	SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, tau,fluxt_zero)
[133]	545	USE icosa
[159]	546	USE omp_para
[133]	547	USE disvert_mod
[159]	548	IMPLICIT NONE
[133]	549	REAL(rstd), INTENT(IN) :: hflux(3iimjjm,llm), wflux(iim*jjm,llm+1)
	550	REAL(rstd), INTENT(INOUT) :: hfluxt(3iimjjm,llm), wfluxt(iim*jjm,llm+1)
	551	REAL(rstd), INTENT(IN) :: tau
[134]	552	LOGICAL, INTENT(INOUT) :: fluxt_zero
[148]	553	INTEGER :: l,i,j,ij
	554
[134]	555	IF(fluxt_zero) THEN
[151]	556
[134]	557	fluxt_zero=.FALSE.
[151]	558
	559	DO l=ll_begin,ll_end
[148]	560	DO j=jj_begin-1,jj_end+1
	561	DO i=ii_begin-1,ii_end+1
	562	ij=(j-1)*iim+i
	563	hfluxt(ij+u_right,l) = tau*hflux(ij+u_right,l)
	564	hfluxt(ij+u_lup,l) = tau*hflux(ij+u_lup,l)
	565	hfluxt(ij+u_ldown,l) = tau*hflux(ij+u_ldown,l)
	566	ENDDO
	567	ENDDO
	568	ENDDO
	569
[159]	570	IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag
	571	DO l=ll_begin,ll_endp1
	572	DO j=jj_begin,jj_end
	573	DO i=ii_begin,ii_end
	574	ij=(j-1)*iim+i
	575	wfluxt(ij,l) = tau*wflux(ij,l)
	576	ENDDO
	577	ENDDO
	578	ENDDO
	579	END IF
[162]	580
[134]	581	ELSE
[151]	582
	583	DO l=ll_begin,ll_end
[148]	584	DO j=jj_begin-1,jj_end+1
	585	DO i=ii_begin-1,ii_end+1
	586	ij=(j-1)*iim+i
	587	hfluxt(ij+u_right,l) = hfluxt(ij+u_right,l)+tau*hflux(ij+u_right,l)
	588	hfluxt(ij+u_lup,l) = hfluxt(ij+u_lup,l)+tau*hflux(ij+u_lup,l)
	589	hfluxt(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)+tau*hflux(ij+u_ldown,l)
	590	ENDDO
	591	ENDDO
	592	ENDDO
	593
[159]	594	IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag
	595	DO l=ll_begin,ll_endp1
	596	DO j=jj_begin,jj_end
	597	DO i=ii_begin,ii_end
	598	ij=(j-1)*iim+i
	599	wfluxt(ij,l) = wfluxt(ij,l)+tau*wflux(ij,l)
	600	ENDDO
	601	ENDDO
	602	ENDDO
	603	END IF
[148]	604
[159]	605	END IF
[151]	606
[133]	607	END SUBROUTINE accumulate_fluxes
[12]	608
[162]	609	FUNCTION maxval_i(p)
	610	USE icosa
	611	IMPLICIT NONE
	612	REAL(rstd), DIMENSION(iim*jjm) :: p
	613	REAL(rstd) :: maxval_i
	614	INTEGER :: j, ij
	615
	616	maxval_i=p((jj_begin-1)*iim+ii_begin)
	617
	618	DO j=jj_begin-1,jj_end+1
	619	ij=(j-1)*iim
	620	maxval_i = MAX(maxval_i, MAXVAL(p(ij+ii_begin:ij+ii_end)))
	621	END DO
	622	END FUNCTION maxval_i
	623
	624	FUNCTION maxval_ik(p)
	625	USE icosa
	626	IMPLICIT NONE
	627	REAL(rstd) :: p(iim*jjm, llm)
	628	REAL(rstd) :: maxval_ik(llm)
	629	INTEGER :: l,j, ij
	630
	631	DO l=1,llm
	632	maxval_ik(l)=p((jj_begin-1)*iim+ii_begin,l)
	633	DO j=jj_begin-1,jj_end+1
	634	ij=(j-1)*iim
	635	maxval_ik(l) = MAX(maxval_ik(l), MAXVAL(p(ij+ii_begin:ij+ii_end,l)))
	636	END DO
	637	END DO
	638	END FUNCTION maxval_ik
	639
[12]	640	END MODULE timeloop_gcm_mod

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

Download in other formats: