Context Navigation

source: codes/icosagcm/trunk/src/caldyn_gcm.f90 @ 171

Last change on this file since 171 was 171, checked in by ymipsl, 11 years ago

XIOS integration -

Compiling with "-with_xios" option. Adapt path to find XIOS library (arch.path)
Retro-compatible with the old output. If xios is not present, dynamico will use the standard writefield function.
Need to have the iodef.xml configuration file in the exec directory

YM

File size: 37.0 KB

Rev	Line
[12]	1	MODULE caldyn_gcm_mod
[19]	2	USE icosa
[151]	3	USE transfert_mod
[132]	4	PRIVATE
	5
	6	INTEGER, PARAMETER :: energy=1, enstrophy=2
[162]	7	TYPE(t_field),POINTER :: f_out_u(:), f_qu(:), f_qv(:)
[159]	8	REAL(rstd),POINTER :: out_u(:,:), p(:,:), qu(:,:)
[17]	9
[50]	10	TYPE(t_field),POINTER :: f_buf_i(:), f_buf_ulon(:), f_buf_ulat(:), f_buf_u3d(:)
	11	TYPE(t_field),POINTER :: f_buf_v(:), f_buf_s(:), f_buf_p(:)
[17]	12
[151]	13	! temporary shared variable for caldyn
	14	TYPE(t_field),POINTER :: f_theta(:)
	15	TYPE(t_field),POINTER :: f_pk(:)
[156]	16	TYPE(t_field),POINTER :: f_geopot(:)
[151]	17	TYPE(t_field),POINTER :: f_wwuu(:)
	18
[159]	19	INTEGER :: caldyn_conserv
[151]	20
[162]	21	TYPE(t_message) :: req_ps, req_mass, req_theta_rhodz, req_u, req_qu
[122]	22
[162]	23	PUBLIC init_caldyn, caldyn_BC, caldyn, write_output_fields, &
[167]	24	req_ps, req_mass
[151]	25
[12]	26	CONTAINS
[15]	27
[98]	28	SUBROUTINE init_caldyn
[50]	29	USE icosa
[122]	30	USE exner_mod
[131]	31	USE mpipara
[50]	32	IMPLICIT NONE
[122]	33	CHARACTER(len=255) :: def
	34
[126]	35	def='enstrophy'
[125]	36	CALL getin('caldyn_conserv',def)
	37	SELECT CASE(TRIM(def))
	38	CASE('energy')
[132]	39	caldyn_conserv=energy
[126]	40	CASE('enstrophy')
[132]	41	caldyn_conserv=enstrophy
[125]	42	CASE DEFAULT
[159]	43	IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_conserv : <', &
	44	TRIM(def),'> options are <energy>, <enstrophy>'
[125]	45	STOP
	46	END SELECT
[134]	47	IF (is_mpi_root) PRINT *, 'caldyn_conserv=',def
[125]	48
[17]	49	CALL allocate_caldyn
[15]	50
	51	END SUBROUTINE init_caldyn
	52
[12]	53	SUBROUTINE allocate_caldyn
[19]	54	USE icosa
[12]	55	IMPLICIT NONE
	56
[151]	57	CALL allocate_field(f_out_u,field_u,type_real,llm)
	58	CALL allocate_field(f_qu,field_u,type_real,llm)
[162]	59	CALL allocate_field(f_qv,field_z,type_real,llm)
[50]	60
[159]	61	CALL allocate_field(f_buf_i, field_t,type_real,llm)
	62	CALL allocate_field(f_buf_p, field_t,type_real,llm+1)
	63	CALL allocate_field(f_buf_u3d, field_t,type_real,3,llm) ! 3D vel at cell centers
[151]	64	CALL allocate_field(f_buf_ulon,field_t,type_real,llm)
	65	CALL allocate_field(f_buf_ulat,field_t,type_real,llm)
[159]	66	CALL allocate_field(f_buf_v, field_z,type_real,llm)
	67	CALL allocate_field(f_buf_s, field_t,type_real)
[50]	68
[159]	69	CALL allocate_field(f_theta, field_t,type_real,llm, name='theta') ! potential temperature
	70	CALL allocate_field(f_pk, field_t,type_real,llm, name='pk')
	71	CALL allocate_field(f_geopot,field_t,type_real,llm+1,name='geopot') ! geopotential
	72	CALL allocate_field(f_wwuu, field_u,type_real,llm+1,name='wwuu')
[151]	73
[12]	74	END SUBROUTINE allocate_caldyn
[157]	75
	76	SUBROUTINE caldyn_BC(f_phis, f_wflux)
	77	USE icosa
	78	USE mpipara
	79	USE omp_para
	80	IMPLICIT NONE
	81	TYPE(t_field),POINTER :: f_phis(:)
	82	TYPE(t_field),POINTER :: f_wflux(:)
	83	REAL(rstd),POINTER :: phis(:)
	84	REAL(rstd),POINTER :: wflux(:,:)
	85	REAL(rstd),POINTER :: geopot(:,:)
	86	REAL(rstd),POINTER :: wwuu(:,:)
	87
	88	INTEGER :: ind,i,j,ij,l
	89
	90	IF (omp_first) THEN
	91	DO ind=1,ndomain
	92	CALL swap_dimensions(ind)
	93	CALL swap_geometry(ind)
	94	geopot=f_geopot(ind)
	95	phis=f_phis(ind)
	96	wflux=f_wflux(ind)
	97	wwuu=f_wwuu(ind)
	98
	99	DO j=jj_begin-1,jj_end+1
	100	DO i=ii_begin-1,ii_end+1
	101	ij=(j-1)*iim+i
	102	! lower BCs : geopot=phis, wflux=0, wwuu=0
	103	geopot(ij,1) = phis(ij)
	104	wflux(ij,1) = 0.
	105	wwuu(ij+u_right,1)=0
	106	wwuu(ij+u_lup,1)=0
	107	wwuu(ij+u_ldown,1)=0
	108	! top BCs : wflux=0, wwuu=0
	109	wflux(ij,llm+1) = 0.
	110	wwuu(ij+u_right,llm+1)=0
	111	wwuu(ij+u_lup,llm+1)=0
	112	wwuu(ij+u_ldown,llm+1)=0
	113	ENDDO
	114	ENDDO
	115	END DO
	116	ENDIF
	117
	118	!$OMP BARRIER
	119	END SUBROUTINE caldyn_BC
[56]	120
[159]	121	SUBROUTINE caldyn(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, &
[162]	122	f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du)
[126]	123	USE icosa
[167]	124	USE disvert_mod, ONLY : caldyn_eta, eta_mass
[126]	125	USE vorticity_mod
	126	USE kinetic_mod
	127	USE theta2theta_rhodz_mod
[131]	128	USE mpipara
[145]	129	USE trace
[151]	130	USE omp_para
[171]	131	USE output_field_mod
[126]	132	IMPLICIT NONE
[129]	133	LOGICAL,INTENT(IN) :: write_out
[126]	134	TYPE(t_field),POINTER :: f_phis(:)
[12]	135	TYPE(t_field),POINTER :: f_ps(:)
[159]	136	TYPE(t_field),POINTER :: f_mass(:)
[126]	137	TYPE(t_field),POINTER :: f_theta_rhodz(:)
	138	TYPE(t_field),POINTER :: f_u(:)
	139	TYPE(t_field),POINTER :: f_q(:)
[134]	140	TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:)
[12]	141	TYPE(t_field),POINTER :: f_dps(:)
[162]	142	TYPE(t_field),POINTER :: f_dmass(:)
[126]	143	TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
	144	TYPE(t_field),POINTER :: f_du(:)
[12]	145
[159]	146	REAL(rstd),POINTER :: ps(:), dps(:)
	147	REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:), dtheta_rhodz(:,:)
[134]	148	REAL(rstd),POINTER :: u(:,:), du(:,:), hflux(:,:), wflux(:,:)
[159]	149	REAL(rstd),POINTER :: qu(:,:)
[162]	150	REAL(rstd),POINTER :: qv(:,:)
[151]	151
	152	! temporary shared variable
	153	REAL(rstd),POINTER :: theta(:,:)
[157]	154	REAL(rstd),POINTER :: pk(:,:)
[156]	155	REAL(rstd),POINTER :: geopot(:,:)
[162]	156	REAL(rstd),POINTER :: convm(:,:)
[151]	157	REAL(rstd),POINTER :: wwuu(:,:)
[157]	158
	159	INTEGER :: ind
[151]	160	LOGICAL,SAVE :: first=.TRUE.
	161	!$OMP THREADPRIVATE(first)
[12]	162
[162]	163	! MPI messages need to be sent at first call to caldyn
	164	! This is needed only once : the next ones will be sent by timeloop
	165	IF (first) THEN
[151]	166	first=.FALSE.
[162]	167	IF(caldyn_eta==eta_mass) THEN
	168	CALL init_message(f_ps,req_i1,req_ps)
	169	ELSE
	170	CALL init_message(f_mass,req_i1,req_mass)
	171	END IF
[151]	172	CALL init_message(f_theta_rhodz,req_i1,req_theta_rhodz)
	173	CALL init_message(f_u,req_e1_vect,req_u)
	174	CALL init_message(f_qu,req_e1_scal,req_qu)
[162]	175	IF(caldyn_eta==eta_mass) THEN
	176	CALL send_message(f_ps,req_ps)
	177	ELSE
	178	CALL send_message(f_mass,req_mass)
	179	END IF
[151]	180	ENDIF
	181
[145]	182	CALL trace_start("caldyn")
[126]	183
[151]	184	CALL send_message(f_u,req_u)
	185	CALL send_message(f_theta_rhodz,req_theta_rhodz)
	186
[126]	187	SELECT CASE(caldyn_conserv)
[132]	188	CASE(energy) ! energy-conserving
[128]	189	DO ind=1,ndomain
	190	CALL swap_dimensions(ind)
	191	CALL swap_geometry(ind)
	192	ps=f_ps(ind)
[157]	193	u=f_u(ind)
	194	theta_rhodz = f_theta_rhodz(ind)
[159]	195	mass=f_mass(ind)
[157]	196	theta = f_theta(ind)
[128]	197	qu=f_qu(ind)
[162]	198	qv=f_qv(ind)
	199	CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu,qv)
[128]	200	ENDDO
	201
[151]	202	CALL send_message(f_qu,req_qu)
[128]	203
	204	DO ind=1,ndomain
	205	CALL swap_dimensions(ind)
	206	CALL swap_geometry(ind)
	207	ps=f_ps(ind)
[157]	208	u=f_u(ind)
[128]	209	theta_rhodz=f_theta_rhodz(ind)
[159]	210	mass=f_mass(ind)
[157]	211	theta = f_theta(ind)
[128]	212	qu=f_qu(ind)
[151]	213	pk = f_pk(ind)
[156]	214	geopot = f_geopot(ind)
[159]	215	CALL compute_geopot(ps,mass,theta, pk,geopot)
[157]	216	hflux=f_hflux(ind)
[162]	217	convm = f_dmass(ind)
[157]	218	dtheta_rhodz=f_dtheta_rhodz(ind)
	219	du=f_du(ind)
[162]	220	CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz,du)
	221	IF(caldyn_eta==eta_mass) THEN
	222	wflux=f_wflux(ind)
	223	wwuu=f_wwuu(ind)
	224	dps=f_dps(ind)
	225	CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
	226	END IF
[128]	227	ENDDO
	228
[132]	229	CASE(enstrophy) ! enstrophy-conserving
[126]	230	DO ind=1,ndomain
	231	CALL swap_dimensions(ind)
	232	CALL swap_geometry(ind)
[157]	233	ps=f_ps(ind)
	234	u=f_u(ind)
	235	theta_rhodz=f_theta_rhodz(ind)
[159]	236	mass=f_mass(ind)
[157]	237	theta = f_theta(ind)
	238	qu=f_qu(ind)
[162]	239	CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu,qv)
[157]	240	pk = f_pk(ind)
[156]	241	geopot = f_geopot(ind)
[159]	242	CALL compute_geopot(ps,mass,theta, pk,geopot)
[134]	243	hflux=f_hflux(ind)
[162]	244	convm = f_dmass(ind)
[126]	245	dtheta_rhodz=f_dtheta_rhodz(ind)
	246	du=f_du(ind)
[162]	247	CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz,du)
	248	IF(caldyn_eta==eta_mass) THEN
	249	wflux=f_wflux(ind)
	250	wwuu=f_wwuu(ind)
	251	dps=f_dps(ind)
	252	CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
	253	END IF
[126]	254	ENDDO
	255
	256	CASE DEFAULT
	257	STOP
	258	END SELECT
[12]	259
[151]	260	!$OMP BARRIER
[129]	261	IF (write_out) THEN
[151]	262
[131]	263	IF (is_mpi_root) PRINT *,'CALL write_output_fields'
[151]	264
	265	! ---> for openMP test to fix later
	266	! CALL write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
	267	! f_buf_i, f_buf_v, f_buf_u3d, f_buf_ulon, f_buf_ulat, f_buf_s, f_buf_p)
[159]	268
[171]	269	CALL output_field("ps",f_ps)
	270	CALL output_field("dps",f_dps)
	271	CALL output_field("mass",f_mass)
	272	CALL output_field("dmass",f_dmass)
	273	CALL output_field("vort",f_qv)
	274	CALL output_field("theta",f_theta)
	275	CALL output_field("exner",f_pk)
	276	CALL output_field("pv",f_qv)
	277
[128]	278	END IF
	279
[126]	280	! CALL check_mass_conservation(f_ps,f_dps)
[145]	281	CALL trace_end("caldyn")
[151]	282	!$OMP BARRIER
[126]	283
	284	END SUBROUTINE caldyn
[128]	285
[162]	286	SUBROUTINE compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu,qv)
[19]	287	USE icosa
[167]	288	USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass
[50]	289	USE exner_mod
[145]	290	USE trace
[151]	291	USE omp_para
[12]	292	IMPLICIT NONE
[128]	293	REAL(rstd),INTENT(IN) :: u(iim3jjm,llm)
	294	REAL(rstd),INTENT(IN) :: ps(iim*jjm)
[157]	295	REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm)
[162]	296	REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm)
[157]	297	REAL(rstd),INTENT(OUT) :: theta(iim*jjm,llm)
[128]	298	REAL(rstd),INTENT(OUT) :: qu(iim3jjm,llm)
[162]	299	REAL(rstd),INTENT(OUT) :: qv(iim2jjm,llm)
[128]	300
	301	INTEGER :: i,j,ij,l
[162]	302	REAL(rstd) :: etav,hv, m
	303	! REAL(rstd) :: qv(2iimjjm,llm) ! potential velocity
[128]	304
[145]	305
[151]	306	CALL trace_start("compute_pvort")
[145]	307
[162]	308	IF(caldyn_eta==eta_mass) THEN
	309	CALL wait_message(req_ps)
	310	ELSE
	311	CALL wait_message(req_mass)
	312	END IF
[157]	313	CALL wait_message(req_theta_rhodz)
[151]	314
[162]	315	IF(caldyn_eta==eta_mass) THEN ! Compute mass & theta
	316	DO l = ll_begin,ll_end
	317	CALL test_message(req_u)
	318	DO j=jj_begin-1,jj_end+1
	319	DO i=ii_begin-1,ii_end+1
	320	ij=(j-1)*iim+i
[167]	321	m = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g
[162]	322	rhodz(ij,l) = m
	323	theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
	324	ENDDO
[128]	325	ENDDO
	326	ENDDO
[162]	327	ELSE ! Compute only theta
	328	DO l = ll_begin,ll_end
	329	CALL test_message(req_u)
	330	DO j=jj_begin-1,jj_end+1
	331	DO i=ii_begin-1,ii_end+1
	332	ij=(j-1)*iim+i
	333	theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
	334	ENDDO
	335	ENDDO
	336	ENDDO
	337	END IF
[151]	338
	339	CALL wait_message(req_u)
[128]	340
[123]	341	!!! Compute shallow-water potential vorticity
[151]	342	DO l = ll_begin,ll_end
	343
[123]	344	DO j=jj_begin-1,jj_end+1
[128]	345	DO i=ii_begin-1,ii_end+1
	346	ij=(j-1)*iim+i
	347
[151]	348	etav= 1./Av(ij+z_up)( ne_rup u(ij+u_rup,l) * de(ij+u_rup) &
	349	+ ne_left * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left) &
	350	- ne_lup * u(ij+u_lup,l) * de(ij+u_lup) )
[123]	351
	352	hv = Riv2(ij,vup) * rhodz(ij,l) &
	353	+ Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l) &
	354	+ Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
	355
	356	qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv
	357
[151]	358	etav = 1./Av(ij+z_down)( ne_ldown u(ij+u_ldown,l) * de(ij+u_ldown) &
	359	+ ne_right * u(ij+t_ldown+u_right,l) * de(ij+t_ldown+u_right) &
	360	- ne_rdown * u(ij+u_rdown,l) * de(ij+u_rdown) )
[123]	361
	362	hv = Riv2(ij,vdown) * rhodz(ij,l) &
	363	+ Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l) &
	364	+ Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
	365
	366	qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv
	367
[12]	368	ENDDO
	369	ENDDO
	370
[126]	371	DO j=jj_begin,jj_end
	372	DO i=ii_begin,ii_end
	373	ij=(j-1)*iim+i
	374	qu(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l))
	375	qu(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l))
	376	qu(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l))
	377	END DO
	378	END DO
	379
	380	ENDDO
	381
[151]	382	CALL trace_end("compute_pvort")
[145]	383
[126]	384	END SUBROUTINE compute_pvort
[125]	385
[159]	386	SUBROUTINE compute_geopot(ps,rhodz,theta,pk,geopot)
[126]	387	USE icosa
	388	USE disvert_mod
	389	USE exner_mod
[145]	390	USE trace
[151]	391	USE omp_para
[126]	392	IMPLICIT NONE
[159]	393	REAL(rstd),INTENT(INOUT) :: ps(iim*jjm)
	394	REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm)
	395	REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm) ! potential temperature
	396	REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) ! Exner function
	397	REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) ! geopotential
[157]	398
	399	INTEGER :: i,j,ij,l
	400	REAL(rstd) :: p_ik, exner_ik
	401
	402	CALL trace_start("compute_geopot")
	403
[159]	404	IF(caldyn_eta==eta_mass) THEN
	405
[157]	406	!!! Compute exner function and geopotential
[159]	407	DO l = 1,llm
	408	!$OMP DO SCHEDULE(STATIC)
	409	DO j=jj_begin-1,jj_end+1
	410	DO i=ii_begin-1,ii_end+1
	411	ij=(j-1)*iim+i
	412	p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later
	413	! p_ik = ptop + g(mass_ak(l)+ mass_bk(l)ps(i,j))
	414	exner_ik = cpp * (p_ik/preff) ** kappa
	415	pk(ij,l) = exner_ik
	416	! specific volume v = kappathetapi/p = dphi/g/rhodz
[162]	417	geopot(ij,l+1) = geopot(ij,l) + (gkappa)rhodz(ij,l)theta(ij,l)exner_ik/p_ik
[159]	418	ENDDO
	419	ENDDO
	420	ENDDO
	421
[167]	422	ELSE
[159]	423	! We are using a Lagrangian vertical coordinate
	424	! Pressure must be computed first top-down (temporarily stored in pk)
	425	! Then Exner pressure and geopotential are computed bottom-up
[167]	426	! Notice that the computation below should work also when caldyn_eta=eta_mass
[159]	427
	428	! uppermost layer
[157]	429	DO j=jj_begin-1,jj_end+1
	430	DO i=ii_begin-1,ii_end+1
	431	ij=(j-1)*iim+i
[159]	432	pk(ij,llm) = ptop + (.5g)rhodz(ij,llm)
	433	END DO
	434	END DO
	435	! other layers
	436	DO l = llm-1, 1, -1
	437	!$OMP DO SCHEDULE(STATIC)
	438	DO j=jj_begin-1,jj_end+1
	439	DO i=ii_begin-1,ii_end+1
	440	ij=(j-1)*iim+i
	441	pk(ij,l) = pk(ij,l+1) + (.5g)(rhodz(ij,l)+rhodz(ij,l+1))
	442	END DO
	443	END DO
	444	END DO
	445	! surface pressure (for diagnostics)
	446	DO j=jj_begin-1,jj_end+1
	447	DO i=ii_begin-1,ii_end+1
	448	ij=(j-1)*iim+i
	449	ps(ij) = pk(ij,1) + (.5g)rhodz(ij,1)
	450	END DO
	451	END DO
[167]	452
	453	IF(boussinesq) THEN ! specific volume 1 = dphi/g/rhodz
	454	DO l = 1,llm
	455	!$OMP DO SCHEDULE(STATIC)
	456	DO j=jj_begin-1,jj_end+1
	457	DO i=ii_begin-1,ii_end+1
	458	ij=(j-1)*iim+i
	459	geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l)
	460	ENDDO
[159]	461	ENDDO
[157]	462	ENDDO
[167]	463	ELSE ! specific volume v = kappathetapi/p = dphi/g/rhodz
	464	DO l = 1,llm
	465	!$OMP DO SCHEDULE(STATIC)
	466	DO j=jj_begin-1,jj_end+1
	467	DO i=ii_begin-1,ii_end+1
	468	ij=(j-1)*iim+i
	469	p_ik = pk(ij,l)
	470	exner_ik = cpp * (p_ik/preff) ** kappa
	471	geopot(ij,l+1) = geopot(ij,l) + (gkappa)rhodz(ij,l)theta(ij,l)exner_ik/p_ik
	472	pk(ij,l) = exner_ik
	473	ENDDO
	474	ENDDO
	475	ENDDO
	476	END IF
[157]	477
[159]	478	END IF
	479
[157]	480	CALL trace_end("compute_geopot")
	481
	482	END SUBROUTINE compute_geopot
	483
[162]	484	SUBROUTINE compute_caldyn_horiz(u,rhodz,qu,theta,pk,geopot, hflux,convm, dtheta_rhodz, du)
[157]	485	USE icosa
	486	USE disvert_mod
	487	USE exner_mod
	488	USE trace
	489	USE omp_para
	490	IMPLICIT NONE
	491	REAL(rstd),INTENT(IN) :: u(iim3jjm,llm)
	492	REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm)
[126]	493	REAL(rstd),INTENT(IN) :: qu(iim3jjm,llm)
[157]	494	REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm) ! potential temperature
[167]	495	REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) ! Exner function
[157]	496	REAL(rstd),INTENT(IN) :: geopot(iim*jjm,llm+1) ! geopotential
[126]	497
[157]	498	REAL(rstd),INTENT(OUT) :: hflux(iim3jjm,llm) ! hflux in kg/s
[162]	499	REAL(rstd),INTENT(OUT) :: convm(iim*jjm,llm) ! mass flux convergence
[126]	500	REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm)
[157]	501	REAL(rstd),INTENT(OUT) :: du(iim3jjm,llm)
[151]	502
	503	REAL(rstd) :: Ftheta(3iimjjm,llm) ! theta flux
[157]	504	REAL(rstd) :: berni(iim*jjm,llm) ! Bernoulli function
[126]	505
[139]	506	INTEGER :: i,j,ij,l
[157]	507	REAL(rstd) :: ww,uu
[139]	508
[157]	509	CALL trace_start("compute_caldyn_horiz")
[126]	510
[151]	511	CALL wait_message(req_theta_rhodz)
	512
	513	DO l = ll_begin, ll_end
[123]	514	!!! Compute mass and theta fluxes
[151]	515	IF (caldyn_conserv==energy) CALL test_message(req_qu)
[12]	516	DO j=jj_begin-1,jj_end+1
	517	DO i=ii_begin-1,ii_end+1
	518	ij=(j-1)*iim+i
[134]	519	hflux(ij+u_right,l)=0.5(rhodz(ij,l)+rhodz(ij+t_right,l))u(ij+u_right,l)*le(ij+u_right)
	520	hflux(ij+u_lup,l)=0.5(rhodz(ij,l)+rhodz(ij+t_lup,l))u(ij+u_lup,l)*le(ij+u_lup)
	521	hflux(ij+u_ldown,l)=0.5(rhodz(ij,l)+rhodz(ij+t_ldown,l))u(ij+u_ldown,l)*le(ij+u_ldown)
[12]	522
[134]	523	Ftheta(ij+u_right,l)=0.5(theta(ij,l)+theta(ij+t_right,l))hflux(ij+u_right,l)
	524	Ftheta(ij+u_lup,l)=0.5(theta(ij,l)+theta(ij+t_lup,l))hflux(ij+u_lup,l)
	525	Ftheta(ij+u_ldown,l)=0.5(theta(ij,l)+theta(ij+t_ldown,l))hflux(ij+u_ldown,l)
[12]	526	ENDDO
	527	ENDDO
[123]	528
	529	!!! compute horizontal divergence of fluxes
[12]	530	DO j=jj_begin,jj_end
	531	DO i=ii_begin,ii_end
	532	ij=(j-1)*iim+i
[162]	533	! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
	534	convm(ij,l)= -1./Ai(ij)(ne_righthflux(ij+u_right,l) + &
[151]	535	ne_rup*hflux(ij+u_rup,l) + &
	536	ne_lup*hflux(ij+u_lup,l) + &
	537	ne_left*hflux(ij+u_left,l) + &
	538	ne_ldown*hflux(ij+u_ldown,l) + &
	539	ne_rdown*hflux(ij+u_rdown,l))
[123]	540
	541	! signe ? attention d (rho theta dz)
[22]	542	! dtheta_rhodz = -div(flux.theta)
[151]	543	dtheta_rhodz(ij,l)=-1./Ai(ij)(ne_rightFtheta(ij+u_right,l) + &
	544	ne_rup*Ftheta(ij+u_rup,l) + &
	545	ne_lup*Ftheta(ij+u_lup,l) + &
	546	ne_left*Ftheta(ij+u_left,l) + &
	547	ne_ldown*Ftheta(ij+u_ldown,l) + &
	548	ne_rdown*Ftheta(ij+u_rdown,l))
[12]	549	ENDDO
	550	ENDDO
[151]	551
[157]	552	END DO
[151]	553
[56]	554	!!! Compute potential vorticity (Coriolis) contribution to du
[157]	555
[128]	556	SELECT CASE(caldyn_conserv)
[132]	557	CASE(energy) ! energy-conserving TRiSK
[12]	558
[151]	559	CALL wait_message(req_qu)
	560
	561	DO l=ll_begin,ll_end
[128]	562	DO j=jj_begin,jj_end
	563	DO i=ii_begin,ii_end
	564	ij=(j-1)*iim+i
[12]	565
[134]	566	uu = wee(ij+u_right,1,1)hflux(ij+u_rup,l)(qu(ij+u_right,l)+qu(ij+u_rup,l))+ &
	567	wee(ij+u_right,2,1)hflux(ij+u_lup,l)(qu(ij+u_right,l)+qu(ij+u_lup,l))+ &
	568	wee(ij+u_right,3,1)hflux(ij+u_left,l)(qu(ij+u_right,l)+qu(ij+u_left,l))+ &
	569	wee(ij+u_right,4,1)hflux(ij+u_ldown,l)(qu(ij+u_right,l)+qu(ij+u_ldown,l))+ &
	570	wee(ij+u_right,5,1)hflux(ij+u_rdown,l)(qu(ij+u_right,l)+qu(ij+u_rdown,l))+ &
	571	wee(ij+u_right,1,2)hflux(ij+t_right+u_ldown,l)(qu(ij+u_right,l)+qu(ij+t_right+u_ldown,l))+ &
	572	wee(ij+u_right,2,2)hflux(ij+t_right+u_rdown,l)(qu(ij+u_right,l)+qu(ij+t_right+u_rdown,l))+ &
	573	wee(ij+u_right,3,2)hflux(ij+t_right+u_right,l)(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))+ &
	574	wee(ij+u_right,4,2)hflux(ij+t_right+u_rup,l)(qu(ij+u_right,l)+qu(ij+t_right+u_rup,l))+ &
	575	wee(ij+u_right,5,2)hflux(ij+t_right+u_lup,l)(qu(ij+u_right,l)+qu(ij+t_right+u_lup,l))
[128]	576	du(ij+u_right,l) = .5*uu/de(ij+u_right)
	577
[134]	578	uu = wee(ij+u_lup,1,1)hflux(ij+u_left,l)(qu(ij+u_lup,l)+qu(ij+u_left,l)) + &
	579	wee(ij+u_lup,2,1)hflux(ij+u_ldown,l)(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) + &
	580	wee(ij+u_lup,3,1)hflux(ij+u_rdown,l)(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) + &
	581	wee(ij+u_lup,4,1)hflux(ij+u_right,l)(qu(ij+u_lup,l)+qu(ij+u_right,l)) + &
	582	wee(ij+u_lup,5,1)hflux(ij+u_rup,l)(qu(ij+u_lup,l)+qu(ij+u_rup,l)) + &
	583	wee(ij+u_lup,1,2)hflux(ij+t_lup+u_right,l)(qu(ij+u_lup,l)+qu(ij+t_lup+u_right,l)) + &
	584	wee(ij+u_lup,2,2)hflux(ij+t_lup+u_rup,l)(qu(ij+u_lup,l)+qu(ij+t_lup+u_rup,l)) + &
	585	wee(ij+u_lup,3,2)hflux(ij+t_lup+u_lup,l)(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) + &
	586	wee(ij+u_lup,4,2)hflux(ij+t_lup+u_left,l)(qu(ij+u_lup,l)+qu(ij+t_lup+u_left,l)) + &
	587	wee(ij+u_lup,5,2)hflux(ij+t_lup+u_ldown,l)(qu(ij+u_lup,l)+qu(ij+t_lup+u_ldown,l))
[128]	588	du(ij+u_lup,l) = .5*uu/de(ij+u_lup)
[12]	589
[128]	590
[134]	591	uu = wee(ij+u_ldown,1,1)hflux(ij+u_rdown,l)(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) + &
	592	wee(ij+u_ldown,2,1)hflux(ij+u_right,l)(qu(ij+u_ldown,l)+qu(ij+u_right,l)) + &
	593	wee(ij+u_ldown,3,1)hflux(ij+u_rup,l)(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) + &
	594	wee(ij+u_ldown,4,1)hflux(ij+u_lup,l)(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) + &
	595	wee(ij+u_ldown,5,1)hflux(ij+u_left,l)(qu(ij+u_ldown,l)+qu(ij+u_left,l)) + &
	596	wee(ij+u_ldown,1,2)hflux(ij+t_ldown+u_lup,l)(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_lup,l)) + &
	597	wee(ij+u_ldown,2,2)hflux(ij+t_ldown+u_left,l)(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_left,l)) + &
	598	wee(ij+u_ldown,3,2)hflux(ij+t_ldown+u_ldown,l)(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) + &
	599	wee(ij+u_ldown,4,2)hflux(ij+t_ldown+u_rdown,l)(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l)) + &
	600	wee(ij+u_ldown,5,2)hflux(ij+t_ldown+u_right,l)(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_right,l))
[128]	601	du(ij+u_ldown,l) = .5*uu/de(ij+u_ldown)
	602
	603	ENDDO
	604	ENDDO
	605	ENDDO
[146]	606
[132]	607	CASE(enstrophy) ! enstrophy-conserving TRiSK
[128]	608
[151]	609	DO l=ll_begin,ll_end
[128]	610	DO j=jj_begin,jj_end
	611	DO i=ii_begin,ii_end
	612	ij=(j-1)*iim+i
[12]	613
[134]	614	uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+ &
	615	wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+ &
	616	wee(ij+u_right,3,1)*hflux(ij+u_left,l)+ &
	617	wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+ &
	618	wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+ &
	619	wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+ &
	620	wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+ &
	621	wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+ &
	622	wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+ &
	623	wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)
[128]	624	du(ij+u_right,l) = qu(ij+u_right,l)*uu/de(ij+u_right)
	625
	626
[134]	627	uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+ &
	628	wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+ &
	629	wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+ &
	630	wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+ &
	631	wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+ &
	632	wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+ &
	633	wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+ &
	634	wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+ &
	635	wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+ &
	636	wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)
[128]	637	du(ij+u_lup,l) = qu(ij+u_lup,l)*uu/de(ij+u_lup)
	638
[134]	639	uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+ &
	640	wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+ &
	641	wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+ &
	642	wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+ &
	643	wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+ &
	644	wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+ &
	645	wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+ &
	646	wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+ &
	647	wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+ &
	648	wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)
[128]	649	du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu/de(ij+u_ldown)
[12]	650
[128]	651	ENDDO
	652	ENDDO
	653	ENDDO
	654
	655	CASE DEFAULT
	656	STOP
	657	END SELECT
[12]	658
	659	!!! Compute bernouilli term = Kinetic Energy + geopotential
[167]	660	IF(boussinesq) THEN
[12]	661
[167]	662	DO l=ll_begin,ll_end
	663	DO j=jj_begin,jj_end
	664	DO i=ii_begin,ii_end
	665	ij=(j-1)*iim+i
	666
	667	berni(ij,l) = pk(ij,l) + &
[12]	668	+ 1/(4Ai(ij))(le(ij+u_right)de(ij+u_right)u(ij+u_right,l)**2 + &
	669	le(ij+u_rup)de(ij+u_rup)u(ij+u_rup,l)**2 + &
	670	le(ij+u_lup)de(ij+u_lup)u(ij+u_lup,l)**2 + &
	671	le(ij+u_left)de(ij+u_left)u(ij+u_left,l)**2 + &
	672	le(ij+u_ldown)de(ij+u_ldown)u(ij+u_ldown,l)**2 + &
	673	le(ij+u_rdown)de(ij+u_rdown)u(ij+u_rdown,l)**2 )
[167]	674	pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
	675	ENDDO
	676	ENDDO
	677	ENDDO
	678
	679	ELSE
	680
	681	DO l=ll_begin,ll_end
	682	DO j=jj_begin,jj_end
	683	DO i=ii_begin,ii_end
	684	ij=(j-1)*iim+i
	685
	686	berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
	687	+ 1/(4Ai(ij))(le(ij+u_right)de(ij+u_right)u(ij+u_right,l)**2 + &
	688	le(ij+u_rup)de(ij+u_rup)u(ij+u_rup,l)**2 + &
	689	le(ij+u_lup)de(ij+u_lup)u(ij+u_lup,l)**2 + &
	690	le(ij+u_left)de(ij+u_left)u(ij+u_left,l)**2 + &
	691	le(ij+u_ldown)de(ij+u_ldown)u(ij+u_ldown,l)**2 + &
	692	le(ij+u_rdown)de(ij+u_rdown)u(ij+u_rdown,l)**2 )
[12]	693
[167]	694	ENDDO
	695	ENDDO
	696	ENDDO
[151]	697
[167]	698	END IF ! Boussinesq/compressible
	699
[157]	700	!!! Add gradients of Bernoulli and Exner functions to du
[167]	701	DO l=ll_begin,ll_end
	702	DO j=jj_begin,jj_end
	703	DO i=ii_begin,ii_end
	704	ij=(j-1)*iim+i
[12]	705
[167]	706	du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right) * ( &
	707	0.5*(theta(ij,l)+theta(ij+t_right,l)) &
	708	( ne_rightpk(ij,l)+ne_left*pk(ij+t_right,l)) &
	709	+ ne_rightberni(ij,l)+ne_leftberni(ij+t_right,l) )
[12]	710
[123]	711
[167]	712	du(ij+u_lup,l) = du(ij+u_lup,l) + 1/de(ij+u_lup) * ( &
	713	0.5*(theta(ij,l)+theta(ij+t_lup,l)) &
	714	( ne_luppk(ij,l)+ne_rdown*pk(ij+t_lup,l)) &
	715	+ ne_lupberni(ij,l)+ne_rdownberni(ij+t_lup,l) )
[123]	716
[167]	717	du(ij+u_ldown,l) = du(ij+u_ldown,l) + 1/de(ij+u_ldown) * ( &
	718	0.5*(theta(ij,l)+theta(ij+t_ldown,l)) &
	719	( ne_ldownpk(ij,l)+ne_rup*pk(ij+t_ldown,l)) &
	720	+ ne_ldownberni(ij,l)+ne_rupberni(ij+t_ldown,l) )
[123]	721
[167]	722	ENDDO
	723	ENDDO
[12]	724	ENDDO
[167]	725
	726	CALL trace_end("compute_caldyn_horiz")
[151]	727
[157]	728	END SUBROUTINE compute_caldyn_horiz
	729
[162]	730	SUBROUTINE compute_caldyn_vert(u,theta,rhodz,convm, wflux,wwuu, dps,dtheta_rhodz,du)
[157]	731	USE icosa
	732	USE disvert_mod
	733	USE exner_mod
	734	USE trace
	735	USE omp_para
	736	IMPLICIT NONE
	737	REAL(rstd),INTENT(IN) :: u(iim3jjm,llm)
	738	REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm)
	739	REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm)
[162]	740	REAL(rstd),INTENT(INOUT) :: convm(iim*jjm,llm) ! mass flux convergence
[157]	741
	742	REAL(rstd),INTENT(OUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s)
	743	REAL(rstd),INTENT(OUT) :: wwuu(iim3jjm,llm+1)
	744	REAL(rstd),INTENT(OUT) :: du(iim3jjm,llm)
	745	REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm)
	746	REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
	747
	748	! temporary variable
	749	INTEGER :: i,j,ij,l
	750	REAL(rstd) :: p_ik, exner_ik
	751
	752	CALL trace_start("compute_caldyn_vert")
	753
	754	!$OMP BARRIER
[162]	755	!!! cumulate mass flux convergence from top to bottom
[157]	756	DO l = llm-1, 1, -1
	757	IF (caldyn_conserv==energy) CALL test_message(req_qu)
	758	!$OMP DO SCHEDULE(STATIC)
[12]	759	DO j=jj_begin,jj_end
	760	DO i=ii_begin,ii_end
[123]	761	ij=(j-1)*iim+i
[162]	762	convm(ij,l) = convm(ij,l) + convm(ij,l+1)
[151]	763	ENDDO
	764	ENDDO
	765	ENDDO
[157]	766
[162]	767	! IMPLICIT FLUSH on convm
[157]	768	!!!!!!!!!!!!!!!!!!!!!!!!!
	769
	770	! compute dps
	771	IF (omp_first) THEN
	772	DO j=jj_begin,jj_end
	773	DO i=ii_begin,ii_end
	774	ij=(j-1)*iim+i
	775	! dps/dt = -int(div flux)dz
[162]	776	dps(ij) = convm(ij,1) * g
[157]	777	ENDDO
	778	ENDDO
	779	ENDIF
[151]	780
[157]	781	!!! Compute vertical mass flux (l=1,llm+1 done by caldyn_BC)
[151]	782	DO l=ll_beginp1,ll_end
[157]	783	IF (caldyn_conserv==energy) CALL test_message(req_qu)
[151]	784	DO j=jj_begin,jj_end
	785	DO i=ii_begin,ii_end
	786	ij=(j-1)*iim+i
[157]	787	! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
	788	! => w>0 for upward transport
[162]	789	wflux( ij, l ) = bp(l) * convm( ij, 1 ) - convm( ij, l )
[151]	790	ENDDO
	791	ENDDO
	792	ENDDO
[157]	793
	794	DO l=ll_begin,ll_endm1
[151]	795	DO j=jj_begin,jj_end
	796	DO i=ii_begin,ii_end
	797	ij=(j-1)*iim+i
[157]	798	dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l ) - 0.5 * ( wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1)))
[151]	799	ENDDO
	800	ENDDO
[157]	801	ENDDO
	802
	803	DO l=ll_beginp1,ll_end
[151]	804	DO j=jj_begin,jj_end
	805	DO i=ii_begin,ii_end
	806	ij=(j-1)*iim+i
[157]	807	dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l ) + 0.5 * ( wflux(ij,l ) * (theta(ij,l-1) + theta(ij,l) ) )
[151]	808	ENDDO
	809	ENDDO
[157]	810	ENDDO
[151]	811
[157]	812	! Compute vertical transport
[151]	813	DO l=ll_beginp1,ll_end
	814	DO j=jj_begin,jj_end
	815	DO i=ii_begin,ii_end
	816	ij=(j-1)*iim+i
	817	wwuu(ij+u_right,l) = 0.5( wflux(ij,l) + wflux(ij+t_right,l)) (u(ij+u_right,l) - u(ij+u_right,l-1))
	818	wwuu(ij+u_lup,l) = 0.5* ( wflux(ij,l) + wflux(ij+t_lup,l)) * (u(ij+u_lup,l) - u(ij+u_lup,l-1))
	819	wwuu(ij+u_ldown,l) = 0.5( wflux(ij,l) + wflux(ij+t_ldown,l)) (u(ij+u_ldown,l) - u(ij+u_ldown,l-1))
	820	ENDDO
	821	ENDDO
	822	ENDDO
[12]	823
[151]	824	!--> flush wwuu
	825	!$OMP BARRIER
[12]	826
[157]	827	! Add vertical transport to du
[151]	828	DO l=ll_begin,ll_end
	829	DO j=jj_begin,jj_end
	830	DO i=ii_begin,ii_end
	831	ij=(j-1)*iim+i
	832	du(ij+u_right, l ) = du(ij+u_right,l) - (wwuu(ij+u_right,l+1)+ wwuu(ij+u_right,l)) / (rhodz(ij,l)+rhodz(ij+t_right,l))
	833	du(ij+u_lup, l ) = du(ij+u_lup,l) - (wwuu(ij+u_lup,l+1) + wwuu(ij+u_lup,l)) / (rhodz(ij,l)+rhodz(ij+t_lup,l))
	834	du(ij+u_ldown, l ) = du(ij+u_ldown,l) - (wwuu(ij+u_ldown,l+1)+ wwuu(ij+u_ldown,l)) / (rhodz(ij,l)+rhodz(ij+t_ldown,l))
[12]	835	ENDDO
	836	ENDDO
	837	ENDDO
[151]	838
[157]	839	CALL trace_end("compute_caldyn_vert")
[145]	840
[157]	841	END SUBROUTINE compute_caldyn_vert
[126]	842
	843	!-------------------------------- Diagnostics ----------------------------
	844
	845	SUBROUTINE check_mass_conservation(f_ps,f_dps)
	846	USE icosa
[131]	847	USE mpipara
[126]	848	IMPLICIT NONE
	849	TYPE(t_field),POINTER :: f_ps(:)
	850	TYPE(t_field),POINTER :: f_dps(:)
	851	REAL(rstd),POINTER :: ps(:)
	852	REAL(rstd),POINTER :: dps(:)
	853	REAL(rstd) :: mass_tot,dmass_tot
	854	INTEGER :: ind,i,j,ij
	855
	856	mass_tot=0
	857	dmass_tot=0
	858
	859	CALL transfert_request(f_dps,req_i1)
	860	CALL transfert_request(f_ps,req_i1)
	861
	862	DO ind=1,ndomain
	863	CALL swap_dimensions(ind)
	864	CALL swap_geometry(ind)
	865
	866	ps=f_ps(ind)
	867	dps=f_dps(ind)
	868
	869	DO j=jj_begin,jj_end
	870	DO i=ii_begin,ii_end
	871	ij=(j-1)*iim+i
	872	IF (domain(ind)%own(i,j)) THEN
	873	mass_tot=mass_tot+ps(ij)*Ai(ij)/g
	874	dmass_tot=dmass_tot+dps(ij)*Ai(ij)/g
	875	ENDIF
	876	ENDDO
	877	ENDDO
	878
	879	ENDDO
[131]	880	IF (is_mpi_root) PRINT*, "mass_tot ", mass_tot," dmass_tot ",dmass_tot
[126]	881
	882	END SUBROUTINE check_mass_conservation
[12]	883
[110]	884	SUBROUTINE write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
[50]	885	f_buf_i, f_buf_v, f_buf_i3, f_buf1_i, f_buf2_i, f_buf_s, f_buf_p)
	886	USE icosa
	887	USE vorticity_mod
	888	USE theta2theta_rhodz_mod
	889	USE pression_mod
[96]	890	USE omega_mod
[50]	891	USE write_field
[97]	892	USE vertical_interp_mod
[151]	893	USE wind_mod
[110]	894	TYPE(t_field),POINTER :: f_ps(:), f_phis(:), f_u(:), f_theta_rhodz(:), f_q(:), f_dps(:), &
[50]	895	f_buf_i(:), f_buf_v(:), f_buf_i3(:), f_buf1_i(:), f_buf2_i(:), f_buf_s(:), f_buf_p(:)
	896
[97]	897	REAL(rstd) :: out_pression_lev
	898	CHARACTER(LEN=255) :: str_pression
[110]	899	CHARACTER(LEN=255) :: physics_type
[97]	900
	901	out_pression_level=0
	902	CALL getin("out_pression_level",out_pression_level)
	903	WRITE(str_pression,*) INT(out_pression_level/100)
	904	str_pression=ADJUSTL(str_pression)
	905
[52]	906	CALL writefield("ps",f_ps)
[151]	907	CALL writefield("dps",f_dps)
	908	CALL writefield("phis",f_phis)
	909	CALL vorticity(f_u,f_buf_v)
	910	CALL writefield("vort",f_buf_v)
[96]	911
[151]	912	CALL w_omega(f_ps, f_u, f_buf_i)
	913	CALL writefield('omega', f_buf_i)
	914	IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
	915	CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
	916	CALL writefield("omega"//TRIM(str_pression),f_buf_s)
	917	ENDIF
[50]	918
	919	! Temperature
[159]	920	! CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ; ! FIXME
	921
[110]	922	CALL getin('physics',physics_type)
	923	IF (TRIM(physics_type)=='dcmip') THEN
	924	CALL Tv2T(f_buf_i,f_q,f_buf1_i)
	925	CALL writefield("T",f_buf1_i)
	926	IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
	927	CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
	928	CALL writefield("T"//TRIM(str_pression),f_buf_s)
	929	ENDIF
	930	ELSE
	931	CALL writefield("T",f_buf_i)
	932	IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
	933	CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
	934	CALL writefield("T"//TRIM(str_pression),f_buf_s)
	935	ENDIF
[97]	936	ENDIF
[110]	937
[50]	938	! velocity components
[151]	939	CALL un2ulonlat(f_u, f_buf1_i, f_buf2_i)
[50]	940	CALL writefield("ulon",f_buf1_i)
	941	CALL writefield("ulat",f_buf2_i)
[97]	942
[104]	943	IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
[97]	944	CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
	945	CALL writefield("ulon"//TRIM(str_pression),f_buf_s)
	946	CALL vertical_interp(f_ps,f_buf2_i,f_buf_s,out_pression_level)
[100]	947	CALL writefield("ulat"//TRIM(str_pression),f_buf_s)
[97]	948	ENDIF
[50]	949
[159]	950	! geopotential ! FIXME
[50]	951	CALL thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_buf_s,f_buf_p,f_buf1_i,f_buf2_i,f_buf_i)
[151]	952	CALL writefield("p",f_buf_p)
[159]	953	CALL writefield("phi",f_geopot) ! geopotential
[151]	954	CALL writefield("theta",f_buf1_i) ! potential temperature
[159]	955	CALL writefield("pk",f_buf2_i) ! Exner pressure
[12]	956
[50]	957	END SUBROUTINE write_output_fields
	958
	959	SUBROUTINE thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_pks,f_p,f_theta,f_pk,f_phi)
	960	USE field_mod
	961	USE pression_mod
	962	USE exner_mod
	963	USE geopotential_mod
	964	USE theta2theta_rhodz_mod
	965	TYPE(t_field), POINTER :: f_ps(:), f_phis(:), f_theta_rhodz(:), & ! IN
	966	f_pks(:), f_p(:), f_theta(:), f_pk(:), f_phi(:) ! OUT
	967	REAL(rstd),POINTER :: pk(:,:), p(:,:), theta(:,:), theta_rhodz(:,:), &
	968	phi(:,:), phis(:), ps(:), pks(:)
	969	INTEGER :: ind
	970
	971	DO ind=1,ndomain
	972	CALL swap_dimensions(ind)
	973	CALL swap_geometry(ind)
	974	ps = f_ps(ind)
	975	p = f_p(ind)
	976	CALL compute_pression(ps,p,0)
	977	pk = f_pk(ind)
	978	pks = f_pks(ind)
	979	CALL compute_exner(ps,p,pks,pk,0)
	980	theta_rhodz = f_theta_rhodz(ind)
	981	theta = f_theta(ind)
	982	CALL compute_theta_rhodz2theta(ps, theta_rhodz,theta,0)
	983	phis = f_phis(ind)
	984	phi = f_phi(ind)
	985	CALL compute_geopotential(phis,pks,pk,theta,phi,0)
	986	END DO
	987
	988	END SUBROUTINE thetarhodz2geopot
	989
[110]	990	SUBROUTINE Tv2T(f_Tv, f_q, f_T)
	991	USE icosa
	992	IMPLICIT NONE
	993	TYPE(t_field), POINTER :: f_TV(:)
	994	TYPE(t_field), POINTER :: f_q(:)
	995	TYPE(t_field), POINTER :: f_T(:)
	996
	997	REAL(rstd),POINTER :: Tv(:,:), q(:,:,:), T(:,:)
	998	INTEGER :: ind
	999
	1000	DO ind=1,ndomain
	1001	CALL swap_dimensions(ind)
	1002	CALL swap_geometry(ind)
	1003	Tv=f_Tv(ind)
	1004	q=f_q(ind)
	1005	T=f_T(ind)
	1006	T=Tv/(1+0.608*q(:,:,1))
	1007	END DO
	1008
	1009	END SUBROUTINE Tv2T
	1010
[12]	1011	END MODULE caldyn_gcm_mod

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

Download in other formats: