Context Navigation

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

Last change on this file since 162 was 162, checked in by dubos, 11 years ago
Lagrangian vertical coordinate tested with test4.1 - 60 MPI procs
File size: 35.5 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, &
	24	req_ps, req_mass, caldyn_eta, eta_mass, eta_lag
[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
[162]	124	USE disvert_mod
[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
[126]	131	IMPLICIT NONE
[129]	132	LOGICAL,INTENT(IN) :: write_out
[126]	133	TYPE(t_field),POINTER :: f_phis(:)
[12]	134	TYPE(t_field),POINTER :: f_ps(:)
[159]	135	TYPE(t_field),POINTER :: f_mass(:)
[126]	136	TYPE(t_field),POINTER :: f_theta_rhodz(:)
	137	TYPE(t_field),POINTER :: f_u(:)
	138	TYPE(t_field),POINTER :: f_q(:)
[134]	139	TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:)
[12]	140	TYPE(t_field),POINTER :: f_dps(:)
[162]	141	TYPE(t_field),POINTER :: f_dmass(:)
[126]	142	TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
	143	TYPE(t_field),POINTER :: f_du(:)
[12]	144
[159]	145	REAL(rstd),POINTER :: ps(:), dps(:)
	146	REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:), dtheta_rhodz(:,:)
[134]	147	REAL(rstd),POINTER :: u(:,:), du(:,:), hflux(:,:), wflux(:,:)
[159]	148	REAL(rstd),POINTER :: qu(:,:)
[162]	149	REAL(rstd),POINTER :: qv(:,:)
[151]	150
	151	! temporary shared variable
	152	REAL(rstd),POINTER :: theta(:,:)
[157]	153	REAL(rstd),POINTER :: pk(:,:)
[156]	154	REAL(rstd),POINTER :: geopot(:,:)
[162]	155	REAL(rstd),POINTER :: convm(:,:)
[151]	156	REAL(rstd),POINTER :: wwuu(:,:)
[157]	157
	158	INTEGER :: ind
[151]	159	LOGICAL,SAVE :: first=.TRUE.
	160	!$OMP THREADPRIVATE(first)
[12]	161
[162]	162	! MPI messages need to be sent at first call to caldyn
	163	! This is needed only once : the next ones will be sent by timeloop
	164	IF (first) THEN
[151]	165	first=.FALSE.
[162]	166	IF(caldyn_eta==eta_mass) THEN
	167	CALL init_message(f_ps,req_i1,req_ps)
	168	ELSE
	169	CALL init_message(f_mass,req_i1,req_mass)
	170	END IF
[151]	171	CALL init_message(f_theta_rhodz,req_i1,req_theta_rhodz)
	172	CALL init_message(f_u,req_e1_vect,req_u)
	173	CALL init_message(f_qu,req_e1_scal,req_qu)
[162]	174	IF(caldyn_eta==eta_mass) THEN
	175	CALL send_message(f_ps,req_ps)
	176	ELSE
	177	CALL send_message(f_mass,req_mass)
	178	END IF
[151]	179	ENDIF
	180
[145]	181	CALL trace_start("caldyn")
[126]	182
[151]	183	CALL send_message(f_u,req_u)
	184	CALL send_message(f_theta_rhodz,req_theta_rhodz)
	185
[126]	186	SELECT CASE(caldyn_conserv)
[132]	187	CASE(energy) ! energy-conserving
[128]	188	DO ind=1,ndomain
	189	CALL swap_dimensions(ind)
	190	CALL swap_geometry(ind)
	191	ps=f_ps(ind)
[157]	192	u=f_u(ind)
	193	theta_rhodz = f_theta_rhodz(ind)
[159]	194	mass=f_mass(ind)
[157]	195	theta = f_theta(ind)
[128]	196	qu=f_qu(ind)
[162]	197	qv=f_qv(ind)
	198	CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu,qv)
[128]	199	ENDDO
	200
[151]	201	CALL send_message(f_qu,req_qu)
[128]	202
	203	DO ind=1,ndomain
	204	CALL swap_dimensions(ind)
	205	CALL swap_geometry(ind)
	206	ps=f_ps(ind)
[157]	207	u=f_u(ind)
[128]	208	theta_rhodz=f_theta_rhodz(ind)
[159]	209	mass=f_mass(ind)
[157]	210	theta = f_theta(ind)
[128]	211	qu=f_qu(ind)
[151]	212	pk = f_pk(ind)
[156]	213	geopot = f_geopot(ind)
[159]	214	CALL compute_geopot(ps,mass,theta, pk,geopot)
[157]	215	hflux=f_hflux(ind)
[162]	216	convm = f_dmass(ind)
[157]	217	dtheta_rhodz=f_dtheta_rhodz(ind)
	218	du=f_du(ind)
[162]	219	CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz,du)
	220	IF(caldyn_eta==eta_mass) THEN
	221	wflux=f_wflux(ind)
	222	wwuu=f_wwuu(ind)
	223	dps=f_dps(ind)
	224	CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
	225	END IF
[128]	226	ENDDO
	227
[132]	228	CASE(enstrophy) ! enstrophy-conserving
[126]	229	DO ind=1,ndomain
	230	CALL swap_dimensions(ind)
	231	CALL swap_geometry(ind)
[157]	232	ps=f_ps(ind)
	233	u=f_u(ind)
	234	theta_rhodz=f_theta_rhodz(ind)
[159]	235	mass=f_mass(ind)
[157]	236	theta = f_theta(ind)
	237	qu=f_qu(ind)
[162]	238	CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu,qv)
[157]	239	pk = f_pk(ind)
[156]	240	geopot = f_geopot(ind)
[159]	241	CALL compute_geopot(ps,mass,theta, pk,geopot)
[134]	242	hflux=f_hflux(ind)
[162]	243	convm = f_dmass(ind)
[126]	244	dtheta_rhodz=f_dtheta_rhodz(ind)
	245	du=f_du(ind)
[162]	246	CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz,du)
	247	IF(caldyn_eta==eta_mass) THEN
	248	wflux=f_wflux(ind)
	249	wwuu=f_wwuu(ind)
	250	dps=f_dps(ind)
	251	CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
	252	END IF
[126]	253	ENDDO
	254
	255	CASE DEFAULT
	256	STOP
	257	END SELECT
[12]	258
[151]	259	!$OMP BARRIER
[129]	260	IF (write_out) THEN
[151]	261
[131]	262	IF (is_mpi_root) PRINT *,'CALL write_output_fields'
[151]	263
	264	! ---> for openMP test to fix later
	265	! CALL write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
	266	! f_buf_i, f_buf_v, f_buf_u3d, f_buf_ulon, f_buf_ulat, f_buf_s, f_buf_p)
[159]	267
[151]	268	CALL writefield("ps",f_ps)
	269	CALL writefield("dps",f_dps)
[159]	270	CALL writefield("mass",f_mass)
[162]	271	CALL writefield("dmass",f_dmass)
[151]	272	CALL writefield("vort",f_qu)
	273	CALL writefield("theta",f_theta)
[162]	274	CALL writefield("exner",f_pk)
	275	CALL writefield("pv",f_qv)
[151]	276
[128]	277	END IF
	278
[126]	279	! CALL check_mass_conservation(f_ps,f_dps)
[145]	280	CALL trace_end("caldyn")
[151]	281	!$OMP BARRIER
[126]	282
	283	END SUBROUTINE caldyn
[128]	284
[162]	285	SUBROUTINE compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu,qv)
[19]	286	USE icosa
[12]	287	USE disvert_mod
[50]	288	USE exner_mod
[145]	289	USE trace
[151]	290	USE omp_para
[12]	291	IMPLICIT NONE
[128]	292	REAL(rstd),INTENT(IN) :: u(iim3jjm,llm)
	293	REAL(rstd),INTENT(IN) :: ps(iim*jjm)
[157]	294	REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm)
[162]	295	REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm)
[157]	296	REAL(rstd),INTENT(OUT) :: theta(iim*jjm,llm)
[128]	297	REAL(rstd),INTENT(OUT) :: qu(iim3jjm,llm)
[162]	298	REAL(rstd),INTENT(OUT) :: qv(iim2jjm,llm)
[128]	299
	300	INTEGER :: i,j,ij,l
[162]	301	REAL(rstd) :: etav,hv, m
	302	! REAL(rstd) :: qv(2iimjjm,llm) ! potential velocity
[128]	303
[145]	304
[151]	305	CALL trace_start("compute_pvort")
[145]	306
[162]	307	IF(caldyn_eta==eta_mass) THEN
	308	CALL wait_message(req_ps)
	309	ELSE
	310	CALL wait_message(req_mass)
	311	END IF
[157]	312	CALL wait_message(req_theta_rhodz)
[151]	313
[162]	314	IF(caldyn_eta==eta_mass) THEN ! Compute mass & theta
	315	DO l = ll_begin,ll_end
	316	CALL test_message(req_u)
	317	DO j=jj_begin-1,jj_end+1
	318	DO i=ii_begin-1,ii_end+1
	319	ij=(j-1)*iim+i
	320	m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g
	321	! m = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g
	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
	422	ELSE
	423
	424	! We are using a Lagrangian vertical coordinate
	425	! Pressure must be computed first top-down (temporarily stored in pk)
	426	! Then Exner pressure and geopotential are computed bottom-up
	427	! Notice that the computation below should work also when caldyn_eta=eta_mass
	428
	429	! uppermost layer
[157]	430	DO j=jj_begin-1,jj_end+1
	431	DO i=ii_begin-1,ii_end+1
	432	ij=(j-1)*iim+i
[159]	433	pk(ij,llm) = ptop + (.5g)rhodz(ij,llm)
	434	END DO
	435	END DO
	436	! other layers
	437	DO l = llm-1, 1, -1
	438	!$OMP DO SCHEDULE(STATIC)
	439	DO j=jj_begin-1,jj_end+1
	440	DO i=ii_begin-1,ii_end+1
	441	ij=(j-1)*iim+i
	442	pk(ij,l) = pk(ij,l+1) + (.5g)(rhodz(ij,l)+rhodz(ij,l+1))
	443	END DO
	444	END DO
	445	END DO
	446	! surface pressure (for diagnostics)
	447	DO j=jj_begin-1,jj_end+1
	448	DO i=ii_begin-1,ii_end+1
	449	ij=(j-1)*iim+i
	450	ps(ij) = pk(ij,1) + (.5g)rhodz(ij,1)
	451	END DO
	452	END DO
	453	DO l = 1,llm
	454	!$OMP DO SCHEDULE(STATIC)
	455	DO j=jj_begin-1,jj_end+1
	456	DO i=ii_begin-1,ii_end+1
	457	ij=(j-1)*iim+i
	458	p_ik = pk(ij,l)
	459	exner_ik = cpp * (p_ik/preff) ** kappa
	460	! specific volume v = kappathetapi/p = dphi/g/rhodz
	461	geopot(ij,l+1) = geopot(ij,l) + (gkappa)rhodz(ij,l)theta(ij,l)exner_ik/p_ik
	462	pk(ij,l) = exner_ik
	463	ENDDO
[157]	464	ENDDO
	465	ENDDO
	466
[159]	467	END IF
	468
[157]	469	CALL trace_end("compute_geopot")
	470
	471	END SUBROUTINE compute_geopot
	472
[162]	473	SUBROUTINE compute_caldyn_horiz(u,rhodz,qu,theta,pk,geopot, hflux,convm, dtheta_rhodz, du)
[157]	474	USE icosa
	475	USE disvert_mod
	476	USE exner_mod
	477	USE trace
	478	USE omp_para
	479	IMPLICIT NONE
	480	REAL(rstd),INTENT(IN) :: u(iim3jjm,llm)
	481	REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm)
[126]	482	REAL(rstd),INTENT(IN) :: qu(iim3jjm,llm)
[157]	483	REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm) ! potential temperature
	484	REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm) ! Exner function
	485	REAL(rstd),INTENT(IN) :: geopot(iim*jjm,llm+1) ! geopotential
[126]	486
[157]	487	REAL(rstd),INTENT(OUT) :: hflux(iim3jjm,llm) ! hflux in kg/s
[162]	488	REAL(rstd),INTENT(OUT) :: convm(iim*jjm,llm) ! mass flux convergence
[126]	489	REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm)
[157]	490	REAL(rstd),INTENT(OUT) :: du(iim3jjm,llm)
[151]	491
	492	REAL(rstd) :: Ftheta(3iimjjm,llm) ! theta flux
[157]	493	REAL(rstd) :: berni(iim*jjm,llm) ! Bernoulli function
[126]	494
[139]	495	INTEGER :: i,j,ij,l
[157]	496	REAL(rstd) :: ww,uu
[139]	497
[157]	498	CALL trace_start("compute_caldyn_horiz")
[126]	499
[151]	500	CALL wait_message(req_theta_rhodz)
	501
	502	DO l = ll_begin, ll_end
[123]	503	!!! Compute mass and theta fluxes
[151]	504	IF (caldyn_conserv==energy) CALL test_message(req_qu)
[12]	505	DO j=jj_begin-1,jj_end+1
	506	DO i=ii_begin-1,ii_end+1
	507	ij=(j-1)*iim+i
[134]	508	hflux(ij+u_right,l)=0.5(rhodz(ij,l)+rhodz(ij+t_right,l))u(ij+u_right,l)*le(ij+u_right)
	509	hflux(ij+u_lup,l)=0.5(rhodz(ij,l)+rhodz(ij+t_lup,l))u(ij+u_lup,l)*le(ij+u_lup)
	510	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]	511
[134]	512	Ftheta(ij+u_right,l)=0.5(theta(ij,l)+theta(ij+t_right,l))hflux(ij+u_right,l)
	513	Ftheta(ij+u_lup,l)=0.5(theta(ij,l)+theta(ij+t_lup,l))hflux(ij+u_lup,l)
	514	Ftheta(ij+u_ldown,l)=0.5(theta(ij,l)+theta(ij+t_ldown,l))hflux(ij+u_ldown,l)
[12]	515	ENDDO
	516	ENDDO
[123]	517
	518	!!! compute horizontal divergence of fluxes
[12]	519	DO j=jj_begin,jj_end
	520	DO i=ii_begin,ii_end
	521	ij=(j-1)*iim+i
[162]	522	! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
	523	convm(ij,l)= -1./Ai(ij)(ne_righthflux(ij+u_right,l) + &
[151]	524	ne_rup*hflux(ij+u_rup,l) + &
	525	ne_lup*hflux(ij+u_lup,l) + &
	526	ne_left*hflux(ij+u_left,l) + &
	527	ne_ldown*hflux(ij+u_ldown,l) + &
	528	ne_rdown*hflux(ij+u_rdown,l))
[123]	529
	530	! signe ? attention d (rho theta dz)
[22]	531	! dtheta_rhodz = -div(flux.theta)
[151]	532	dtheta_rhodz(ij,l)=-1./Ai(ij)(ne_rightFtheta(ij+u_right,l) + &
	533	ne_rup*Ftheta(ij+u_rup,l) + &
	534	ne_lup*Ftheta(ij+u_lup,l) + &
	535	ne_left*Ftheta(ij+u_left,l) + &
	536	ne_ldown*Ftheta(ij+u_ldown,l) + &
	537	ne_rdown*Ftheta(ij+u_rdown,l))
[12]	538	ENDDO
	539	ENDDO
[151]	540
[157]	541	END DO
[151]	542
[56]	543	!!! Compute potential vorticity (Coriolis) contribution to du
[157]	544
[128]	545	SELECT CASE(caldyn_conserv)
[132]	546	CASE(energy) ! energy-conserving TRiSK
[12]	547
[151]	548	CALL wait_message(req_qu)
	549
	550	DO l=ll_begin,ll_end
[128]	551	DO j=jj_begin,jj_end
	552	DO i=ii_begin,ii_end
	553	ij=(j-1)*iim+i
[12]	554
[134]	555	uu = wee(ij+u_right,1,1)hflux(ij+u_rup,l)(qu(ij+u_right,l)+qu(ij+u_rup,l))+ &
	556	wee(ij+u_right,2,1)hflux(ij+u_lup,l)(qu(ij+u_right,l)+qu(ij+u_lup,l))+ &
	557	wee(ij+u_right,3,1)hflux(ij+u_left,l)(qu(ij+u_right,l)+qu(ij+u_left,l))+ &
	558	wee(ij+u_right,4,1)hflux(ij+u_ldown,l)(qu(ij+u_right,l)+qu(ij+u_ldown,l))+ &
	559	wee(ij+u_right,5,1)hflux(ij+u_rdown,l)(qu(ij+u_right,l)+qu(ij+u_rdown,l))+ &
	560	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))+ &
	561	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))+ &
	562	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))+ &
	563	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))+ &
	564	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]	565	du(ij+u_right,l) = .5*uu/de(ij+u_right)
	566
[134]	567	uu = wee(ij+u_lup,1,1)hflux(ij+u_left,l)(qu(ij+u_lup,l)+qu(ij+u_left,l)) + &
	568	wee(ij+u_lup,2,1)hflux(ij+u_ldown,l)(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) + &
	569	wee(ij+u_lup,3,1)hflux(ij+u_rdown,l)(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) + &
	570	wee(ij+u_lup,4,1)hflux(ij+u_right,l)(qu(ij+u_lup,l)+qu(ij+u_right,l)) + &
	571	wee(ij+u_lup,5,1)hflux(ij+u_rup,l)(qu(ij+u_lup,l)+qu(ij+u_rup,l)) + &
	572	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)) + &
	573	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)) + &
	574	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)) + &
	575	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)) + &
	576	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]	577	du(ij+u_lup,l) = .5*uu/de(ij+u_lup)
[12]	578
[128]	579
[134]	580	uu = wee(ij+u_ldown,1,1)hflux(ij+u_rdown,l)(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) + &
	581	wee(ij+u_ldown,2,1)hflux(ij+u_right,l)(qu(ij+u_ldown,l)+qu(ij+u_right,l)) + &
	582	wee(ij+u_ldown,3,1)hflux(ij+u_rup,l)(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) + &
	583	wee(ij+u_ldown,4,1)hflux(ij+u_lup,l)(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) + &
	584	wee(ij+u_ldown,5,1)hflux(ij+u_left,l)(qu(ij+u_ldown,l)+qu(ij+u_left,l)) + &
	585	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)) + &
	586	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)) + &
	587	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)) + &
	588	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)) + &
	589	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]	590	du(ij+u_ldown,l) = .5*uu/de(ij+u_ldown)
	591
	592	ENDDO
	593	ENDDO
	594	ENDDO
[146]	595
[132]	596	CASE(enstrophy) ! enstrophy-conserving TRiSK
[128]	597
[151]	598	DO l=ll_begin,ll_end
[128]	599	DO j=jj_begin,jj_end
	600	DO i=ii_begin,ii_end
	601	ij=(j-1)*iim+i
[12]	602
[134]	603	uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+ &
	604	wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+ &
	605	wee(ij+u_right,3,1)*hflux(ij+u_left,l)+ &
	606	wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+ &
	607	wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+ &
	608	wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+ &
	609	wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+ &
	610	wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+ &
	611	wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+ &
	612	wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)
[128]	613	du(ij+u_right,l) = qu(ij+u_right,l)*uu/de(ij+u_right)
	614
	615
[134]	616	uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+ &
	617	wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+ &
	618	wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+ &
	619	wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+ &
	620	wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+ &
	621	wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+ &
	622	wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+ &
	623	wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+ &
	624	wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+ &
	625	wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)
[128]	626	du(ij+u_lup,l) = qu(ij+u_lup,l)*uu/de(ij+u_lup)
	627
[134]	628	uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+ &
	629	wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+ &
	630	wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+ &
	631	wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+ &
	632	wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+ &
	633	wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+ &
	634	wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+ &
	635	wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+ &
	636	wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+ &
	637	wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)
[128]	638	du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu/de(ij+u_ldown)
[12]	639
[128]	640	ENDDO
	641	ENDDO
	642	ENDDO
	643
	644	CASE DEFAULT
	645	STOP
	646	END SELECT
[12]	647
	648	!!! Compute bernouilli term = Kinetic Energy + geopotential
[151]	649	DO l=ll_begin,ll_end
[12]	650	DO j=jj_begin,jj_end
	651	DO i=ii_begin,ii_end
	652	ij=(j-1)*iim+i
	653
[156]	654	berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
[12]	655	+ 1/(4Ai(ij))(le(ij+u_right)de(ij+u_right)u(ij+u_right,l)**2 + &
	656	le(ij+u_rup)de(ij+u_rup)u(ij+u_rup,l)**2 + &
	657	le(ij+u_lup)de(ij+u_lup)u(ij+u_lup,l)**2 + &
	658	le(ij+u_left)de(ij+u_left)u(ij+u_left,l)**2 + &
	659	le(ij+u_ldown)de(ij+u_ldown)u(ij+u_ldown,l)**2 + &
	660	le(ij+u_rdown)de(ij+u_rdown)u(ij+u_rdown,l)**2 )
	661
	662	ENDDO
	663	ENDDO
	664	ENDDO
[151]	665
[12]	666
[157]	667	!!! Add gradients of Bernoulli and Exner functions to du
[151]	668
	669	DO l=ll_begin,ll_end
[12]	670	DO j=jj_begin,jj_end
	671	DO i=ii_begin,ii_end
	672	ij=(j-1)*iim+i
	673
[151]	674	du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right) * ( &
	675	0.5*(theta(ij,l)+theta(ij+t_right,l)) &
	676	( ne_rightpk(ij,l)+ne_left*pk(ij+t_right,l)) &
	677	+ ne_rightberni(ij,l)+ne_leftberni(ij+t_right,l) )
[12]	678
[123]	679
[151]	680	du(ij+u_lup,l) = du(ij+u_lup,l) + 1/de(ij+u_lup) * ( &
	681	0.5*(theta(ij,l)+theta(ij+t_lup,l)) &
	682	( ne_luppk(ij,l)+ne_rdown*pk(ij+t_lup,l)) &
	683	+ ne_lupberni(ij,l)+ne_rdownberni(ij+t_lup,l) )
[123]	684
[151]	685	du(ij+u_ldown,l) = du(ij+u_ldown,l) + 1/de(ij+u_ldown) * ( &
	686	0.5*(theta(ij,l)+theta(ij+t_ldown,l)) &
	687	( ne_ldownpk(ij,l)+ne_rup*pk(ij+t_ldown,l)) &
	688	+ ne_ldownberni(ij,l)+ne_rupberni(ij+t_ldown,l) )
[123]	689
[12]	690	ENDDO
	691	ENDDO
	692	ENDDO
[151]	693
[157]	694	CALL trace_end("compute_caldyn_horiz")
	695
	696	END SUBROUTINE compute_caldyn_horiz
	697
[162]	698	SUBROUTINE compute_caldyn_vert(u,theta,rhodz,convm, wflux,wwuu, dps,dtheta_rhodz,du)
[157]	699	USE icosa
	700	USE disvert_mod
	701	USE exner_mod
	702	USE trace
	703	USE omp_para
	704	IMPLICIT NONE
	705	REAL(rstd),INTENT(IN) :: u(iim3jjm,llm)
	706	REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm)
	707	REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm)
[162]	708	REAL(rstd),INTENT(INOUT) :: convm(iim*jjm,llm) ! mass flux convergence
[157]	709
	710	REAL(rstd),INTENT(OUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s)
	711	REAL(rstd),INTENT(OUT) :: wwuu(iim3jjm,llm+1)
	712	REAL(rstd),INTENT(OUT) :: du(iim3jjm,llm)
	713	REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm)
	714	REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
	715
	716	! temporary variable
	717	INTEGER :: i,j,ij,l
	718	REAL(rstd) :: p_ik, exner_ik
	719
	720	CALL trace_start("compute_caldyn_vert")
	721
	722	!$OMP BARRIER
[162]	723	!!! cumulate mass flux convergence from top to bottom
[157]	724	DO l = llm-1, 1, -1
	725	IF (caldyn_conserv==energy) CALL test_message(req_qu)
	726	!$OMP DO SCHEDULE(STATIC)
[12]	727	DO j=jj_begin,jj_end
	728	DO i=ii_begin,ii_end
[123]	729	ij=(j-1)*iim+i
[162]	730	convm(ij,l) = convm(ij,l) + convm(ij,l+1)
[151]	731	ENDDO
	732	ENDDO
	733	ENDDO
[157]	734
[162]	735	! IMPLICIT FLUSH on convm
[157]	736	!!!!!!!!!!!!!!!!!!!!!!!!!
	737
	738	! compute dps
	739	IF (omp_first) THEN
	740	DO j=jj_begin,jj_end
	741	DO i=ii_begin,ii_end
	742	ij=(j-1)*iim+i
	743	! dps/dt = -int(div flux)dz
[162]	744	dps(ij) = convm(ij,1) * g
[157]	745	ENDDO
	746	ENDDO
	747	ENDIF
[151]	748
[157]	749	!!! Compute vertical mass flux (l=1,llm+1 done by caldyn_BC)
[151]	750	DO l=ll_beginp1,ll_end
[157]	751	IF (caldyn_conserv==energy) CALL test_message(req_qu)
[151]	752	DO j=jj_begin,jj_end
	753	DO i=ii_begin,ii_end
	754	ij=(j-1)*iim+i
[157]	755	! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
	756	! => w>0 for upward transport
[162]	757	wflux( ij, l ) = bp(l) * convm( ij, 1 ) - convm( ij, l )
[151]	758	ENDDO
	759	ENDDO
	760	ENDDO
[157]	761
	762	DO l=ll_begin,ll_endm1
[151]	763	DO j=jj_begin,jj_end
	764	DO i=ii_begin,ii_end
	765	ij=(j-1)*iim+i
[157]	766	dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l ) - 0.5 * ( wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1)))
[151]	767	ENDDO
	768	ENDDO
[157]	769	ENDDO
	770
	771	DO l=ll_beginp1,ll_end
[151]	772	DO j=jj_begin,jj_end
	773	DO i=ii_begin,ii_end
	774	ij=(j-1)*iim+i
[157]	775	dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l ) + 0.5 * ( wflux(ij,l ) * (theta(ij,l-1) + theta(ij,l) ) )
[151]	776	ENDDO
	777	ENDDO
[157]	778	ENDDO
[151]	779
[157]	780	! Compute vertical transport
[151]	781	DO l=ll_beginp1,ll_end
	782	DO j=jj_begin,jj_end
	783	DO i=ii_begin,ii_end
	784	ij=(j-1)*iim+i
	785	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))
	786	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))
	787	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))
	788	ENDDO
	789	ENDDO
	790	ENDDO
[12]	791
[151]	792	!--> flush wwuu
	793	!$OMP BARRIER
[12]	794
[157]	795	! Add vertical transport to du
[151]	796	DO l=ll_begin,ll_end
	797	DO j=jj_begin,jj_end
	798	DO i=ii_begin,ii_end
	799	ij=(j-1)*iim+i
	800	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))
	801	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))
	802	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]	803	ENDDO
	804	ENDDO
	805	ENDDO
[151]	806
[157]	807	CALL trace_end("compute_caldyn_vert")
[145]	808
[157]	809	END SUBROUTINE compute_caldyn_vert
[126]	810
	811	!-------------------------------- Diagnostics ----------------------------
	812
	813	SUBROUTINE check_mass_conservation(f_ps,f_dps)
	814	USE icosa
[131]	815	USE mpipara
[126]	816	IMPLICIT NONE
	817	TYPE(t_field),POINTER :: f_ps(:)
	818	TYPE(t_field),POINTER :: f_dps(:)
	819	REAL(rstd),POINTER :: ps(:)
	820	REAL(rstd),POINTER :: dps(:)
	821	REAL(rstd) :: mass_tot,dmass_tot
	822	INTEGER :: ind,i,j,ij
	823
	824	mass_tot=0
	825	dmass_tot=0
	826
	827	CALL transfert_request(f_dps,req_i1)
	828	CALL transfert_request(f_ps,req_i1)
	829
	830	DO ind=1,ndomain
	831	CALL swap_dimensions(ind)
	832	CALL swap_geometry(ind)
	833
	834	ps=f_ps(ind)
	835	dps=f_dps(ind)
	836
	837	DO j=jj_begin,jj_end
	838	DO i=ii_begin,ii_end
	839	ij=(j-1)*iim+i
	840	IF (domain(ind)%own(i,j)) THEN
	841	mass_tot=mass_tot+ps(ij)*Ai(ij)/g
	842	dmass_tot=dmass_tot+dps(ij)*Ai(ij)/g
	843	ENDIF
	844	ENDDO
	845	ENDDO
	846
	847	ENDDO
[131]	848	IF (is_mpi_root) PRINT*, "mass_tot ", mass_tot," dmass_tot ",dmass_tot
[126]	849
	850	END SUBROUTINE check_mass_conservation
[12]	851
[110]	852	SUBROUTINE write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
[50]	853	f_buf_i, f_buf_v, f_buf_i3, f_buf1_i, f_buf2_i, f_buf_s, f_buf_p)
	854	USE icosa
	855	USE vorticity_mod
	856	USE theta2theta_rhodz_mod
	857	USE pression_mod
[96]	858	USE omega_mod
[50]	859	USE write_field
[97]	860	USE vertical_interp_mod
[151]	861	USE wind_mod
[110]	862	TYPE(t_field),POINTER :: f_ps(:), f_phis(:), f_u(:), f_theta_rhodz(:), f_q(:), f_dps(:), &
[50]	863	f_buf_i(:), f_buf_v(:), f_buf_i3(:), f_buf1_i(:), f_buf2_i(:), f_buf_s(:), f_buf_p(:)
	864
[97]	865	REAL(rstd) :: out_pression_lev
	866	CHARACTER(LEN=255) :: str_pression
[110]	867	CHARACTER(LEN=255) :: physics_type
[97]	868
	869	out_pression_level=0
	870	CALL getin("out_pression_level",out_pression_level)
	871	WRITE(str_pression,*) INT(out_pression_level/100)
	872	str_pression=ADJUSTL(str_pression)
	873
[52]	874	CALL writefield("ps",f_ps)
[151]	875	CALL writefield("dps",f_dps)
	876	CALL writefield("phis",f_phis)
	877	CALL vorticity(f_u,f_buf_v)
	878	CALL writefield("vort",f_buf_v)
[96]	879
[151]	880	CALL w_omega(f_ps, f_u, f_buf_i)
	881	CALL writefield('omega', f_buf_i)
	882	IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
	883	CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
	884	CALL writefield("omega"//TRIM(str_pression),f_buf_s)
	885	ENDIF
[50]	886
	887	! Temperature
[159]	888	! CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ; ! FIXME
	889
[110]	890	CALL getin('physics',physics_type)
	891	IF (TRIM(physics_type)=='dcmip') THEN
	892	CALL Tv2T(f_buf_i,f_q,f_buf1_i)
	893	CALL writefield("T",f_buf1_i)
	894	IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
	895	CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
	896	CALL writefield("T"//TRIM(str_pression),f_buf_s)
	897	ENDIF
	898	ELSE
	899	CALL writefield("T",f_buf_i)
	900	IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
	901	CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
	902	CALL writefield("T"//TRIM(str_pression),f_buf_s)
	903	ENDIF
[97]	904	ENDIF
[110]	905
[50]	906	! velocity components
[151]	907	CALL un2ulonlat(f_u, f_buf1_i, f_buf2_i)
[50]	908	CALL writefield("ulon",f_buf1_i)
	909	CALL writefield("ulat",f_buf2_i)
[97]	910
[104]	911	IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
[97]	912	CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
	913	CALL writefield("ulon"//TRIM(str_pression),f_buf_s)
	914	CALL vertical_interp(f_ps,f_buf2_i,f_buf_s,out_pression_level)
[100]	915	CALL writefield("ulat"//TRIM(str_pression),f_buf_s)
[97]	916	ENDIF
[50]	917
[159]	918	! geopotential ! FIXME
[50]	919	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]	920	CALL writefield("p",f_buf_p)
[159]	921	CALL writefield("phi",f_geopot) ! geopotential
[151]	922	CALL writefield("theta",f_buf1_i) ! potential temperature
[159]	923	CALL writefield("pk",f_buf2_i) ! Exner pressure
[12]	924
[50]	925	END SUBROUTINE write_output_fields
	926
	927	SUBROUTINE thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_pks,f_p,f_theta,f_pk,f_phi)
	928	USE field_mod
	929	USE pression_mod
	930	USE exner_mod
	931	USE geopotential_mod
	932	USE theta2theta_rhodz_mod
	933	TYPE(t_field), POINTER :: f_ps(:), f_phis(:), f_theta_rhodz(:), & ! IN
	934	f_pks(:), f_p(:), f_theta(:), f_pk(:), f_phi(:) ! OUT
	935	REAL(rstd),POINTER :: pk(:,:), p(:,:), theta(:,:), theta_rhodz(:,:), &
	936	phi(:,:), phis(:), ps(:), pks(:)
	937	INTEGER :: ind
	938
	939	DO ind=1,ndomain
	940	CALL swap_dimensions(ind)
	941	CALL swap_geometry(ind)
	942	ps = f_ps(ind)
	943	p = f_p(ind)
	944	CALL compute_pression(ps,p,0)
	945	pk = f_pk(ind)
	946	pks = f_pks(ind)
	947	CALL compute_exner(ps,p,pks,pk,0)
	948	theta_rhodz = f_theta_rhodz(ind)
	949	theta = f_theta(ind)
	950	CALL compute_theta_rhodz2theta(ps, theta_rhodz,theta,0)
	951	phis = f_phis(ind)
	952	phi = f_phi(ind)
	953	CALL compute_geopotential(phis,pks,pk,theta,phi,0)
	954	END DO
	955
	956	END SUBROUTINE thetarhodz2geopot
	957
[110]	958	SUBROUTINE Tv2T(f_Tv, f_q, f_T)
	959	USE icosa
	960	IMPLICIT NONE
	961	TYPE(t_field), POINTER :: f_TV(:)
	962	TYPE(t_field), POINTER :: f_q(:)
	963	TYPE(t_field), POINTER :: f_T(:)
	964
	965	REAL(rstd),POINTER :: Tv(:,:), q(:,:,:), T(:,:)
	966	INTEGER :: ind
	967
	968	DO ind=1,ndomain
	969	CALL swap_dimensions(ind)
	970	CALL swap_geometry(ind)
	971	Tv=f_Tv(ind)
	972	q=f_q(ind)
	973	T=f_T(ind)
	974	T=Tv/(1+0.608*q(:,:,1))
	975	END DO
	976
	977	END SUBROUTINE Tv2T
	978
[12]	979	END MODULE caldyn_gcm_mod

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

Download in other formats: