Context Navigation

source: codes/icosagcm/devel/src/diagnostics/check_conserve.f90 @ 609

Last change on this file since 609 was 609, checked in by dubos, 7 years ago
devel : backported commits 606 and 607 from trunk
File size: 14.3 KB

Rev	Line
[149]	1	MODULE check_conserve_mod
[151]	2	USE icosa
	3	IMPLICIT NONE
[149]	4
[365]	5	PRIVATE
	6
	7	INTEGER, PARAMETER, PUBLIC :: check_basic=1, check_detailed=2, &
	8	AAM_dissip=1, AAM_dyn=2, AAM_phys=3
	9	INTEGER, SAVE :: check_type
[326]	10
[365]	11	TYPE(t_field),POINTER,SAVE :: f_pk(:), f_pks(:), f_p(:)
	12	TYPE(t_field),POINTER,SAVE :: f_vort(:)
	13	TYPE(t_field),POINTER,SAVE :: f_rhodz(:)
	14
	15	REAL(rstd),SAVE :: AAM_mass, AAM_mass_old, AAM_vel, AAM_vel_old, &
	16	AAM_mass_source(3), AAM_vel_source(3) ! read/written only IF is_master
	17	REAL(rstd),SAVE :: mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0
	18	!$OMP THREADPRIVATE(check_type, mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0)
	19
	20	PUBLIC :: init_check_conserve, check_conserve_detailed, check_conserve
[151]	21	CONTAINS
[149]	22
	23	!---------------------------------------------------------------------
	24
[151]	25	SUBROUTINE init_check_conserve
[326]	26	USE getin_mod
	27	USE omp_para, ONLY : is_master
	28	CHARACTER(LEN=255) :: check_type_str
[149]	29	CALL allocate_field(f_pk,field_t,type_real,llm)
	30	CALL allocate_field(f_p,field_t,type_real,llm+1)
	31	CALL allocate_field(f_pks,field_t,type_real)
	32	CALL allocate_field(f_rhodz,field_t,type_real,llm)
	33	CALL allocate_field(f_vort,field_z,type_real,llm)
[326]	34
	35	check_type_str='basic'
	36	CALL getin("check_conservation",check_type_str)
	37	SELECT CASE(TRIM(check_type_str))
	38	CASE('basic')
	39	check_type=check_basic
	40	CASE('detailed')
	41	check_type=check_detailed
	42	CASE DEFAULT
	43	IF(is_master) PRINT*,'Bad selector <', TRIM(check_type_str), &
	44	'> for variable check_conservation : options are <basic>, <detailed>'
	45	STOP
	46	END SELECT
[151]	47	END SUBROUTINE init_check_conserve
[149]	48
[326]	49	!--------------------------------- Basic check --------------------------------
	50
[151]	51	SUBROUTINE check_conserve(f_ps,f_dps,f_ue,f_theta_rhodz,f_phis,it)
	52	USE pression_mod
	53	USE vorticity_mod
	54	USE caldyn_gcm_mod
	55	USE exner_mod
[365]	56	! USE mpipara, ONLY : is_mpi_master, comm_icosa
	57	USE omp_para, ONLY : is_master
[151]	58	TYPE(t_field),POINTER :: f_ps(:)
	59	TYPE(t_field),POINTER :: f_dps(:)
	60	TYPE(t_field),POINTER :: f_ue(:)
	61	TYPE(t_field),POINTER :: f_theta_rhodz(:)
	62	TYPE(t_field),POINTER :: f_phis(:)
[326]	63	INTEGER, INTENT(IN) :: it
[160]	64
[151]	65	REAL(rstd),POINTER :: p(:,:),rhodz(:,:)
[365]	66	INTEGER :: ind,ierr
[326]	67	REAL(rstd) :: mtot, angtot, rmsdpdt
[365]	68	REAL(rstd) :: etot, stot, ang_mass, ang_vel, rmsvtot, ztot
[151]	69
[295]	70	CALL transfert_request(f_ue,req_e1_vect)
[151]	71	CALL pression(f_ps,f_p)
	72
	73	DO ind=1,ndomain
[186]	74	IF (.NOT. assigned_domain(ind)) CYCLE
[151]	75	CALL swap_dimensions(ind)
	76	CALL swap_geometry(ind)
[160]	77	p=f_p(ind)
	78	rhodz=f_rhodz(ind)
[151]	79	CALL compute_rhodz(p,rhodz)
[160]	80	END DO
[151]	81
	82	CALL vorticity(f_ue,f_vort)
[160]	83	CALL check_mass_conserve(f_ps,f_dps,mtot,rmsdpdt)
[295]	84	CALL check_PV(ztot)
[151]	85	CALL exner(f_ps,f_p,f_pks,f_pk)
[365]	86	CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vel, rmsvtot)
[160]	87
[365]	88	IF (is_master) THEN ! is_master = is_omp_master && is_mpi_master
	89	! AAM_mass and AAM_vel must be read/written only IF is_master
	90	AAM_mass = ang_mass
	91	AAM_vel = ang_vel
	92	angtot = ang_mass + ang_vel
[295]	93	IF ( it == itau0 ) THEN
[160]	94	ztot0 = ztot
	95	mtot0 = mtot
	96	etot0 = etot
[326]	97	angtot0 = angtot
[160]	98	stot0 = stot
[365]	99	AAM_mass_old=AAM_mass
	100	AAM_vel_old=AAM_vel
	101	AAM_mass_source=0.
	102	AAM_vel_source=0.
[160]	103	END IF
[295]	104	rmsvtot=SQRT(rmsvtot/mtot)
[160]	105	ztot=ztot/ztot0-1. ; mtot=mtot/mtot0-1.
[295]	106	etot=etot/etot0-1. ; angtot=angtot/angtot0-1. ; stot=stot/stot0-1.
[160]	107	rmsdpdt= daysec1.e-2sqrt(rmsdpdt/ncell_glo)
	108
[151]	109	OPEN(134,file="checkconsicosa.txt",position='append')
[295]	110	WRITE(134,4000)mtot,rmsdpdt,etot,ztot,stot,rmsvtot,angtot
	111	WRITE(134,*)mtot,rmsdpdt,etot,ztot,stot,rmsvtot,angtot
[151]	112	WRITE(134,*)"=================================================="
[295]	113	WRITE(*,4000)mtot,rmsdpdt,etot,ztot,stot,rmsvtot,angtot
[160]	114
[186]	115	4000 FORMAT(10x,'masse',5x,'rmsdpdt',5x,'energie',5x,'enstrophie' &
[365]	116	,5x,'entropie',5x,'rmsv',5x,'mt.ang',/,'GLOB ' &
[295]	117	,e10.3,e13.6,5e13.3/)
[365]	118	CLOSE(134)
	119
	120	IF(check_type == check_detailed) THEN
	121	!ang_mass = AAM_mass-AAM_mass_old
	122	!WRITE(*,'(A,3E12.4)') 'AAM_mass sanity check', SUM(AAM_mass_source), ang_mass, SUM(AAM_mass_source)-ang_mass
	123	!ang_vel = AAM_vel-AAM_vel_old
	124	!WRITE(*,'(A,3E12.4)') 'AAM_vel sanity check', SUM(AAM_vel_source), ang_vel, SUM(AAM_vel_source)-ang_vel
	125	WRITE(,'(A,6E12.4)') 'AAM_mass : time,old,new,dissip,dyn,phys', dtit, AAM_mass_old, AAM_mass, &
	126	AAM_mass_source(AAM_dissip), AAM_mass_source(AAM_dyn), AAM_mass_source(AAM_phys)
	127	WRITE(,'(A,6E12.4)') 'AAM_vel : time,old,new,dissip,dyn,phys', dtit, AAM_vel_old, AAM_vel, &
	128	AAM_vel_source(AAM_dissip), AAM_vel_source(AAM_dyn), AAM_vel_source(AAM_phys)
[609]	129	WRITE(,'(A,6E12.4)') 'AAM_tot : time,old,new,dissip,dyn,phys', dtit, &
	130	AAM_vel_old+AAM_mass_old, AAM_vel+AAM_mass, &
	131	AAM_vel_source(AAM_dissip)+AAM_mass_source(AAM_dissip), &
	132	AAM_vel_source(AAM_dyn) +AAM_mass_source(AAM_dyn), &
	133	AAM_vel_source(AAM_phys) +AAM_mass_source(AAM_phys)
[365]	134	AAM_mass_old=AAM_mass
	135	AAM_vel_old=AAM_vel
	136	AAM_mass_source=0.
	137	AAM_vel_source=0.
	138	END IF
[160]	139	END IF
[365]	140
[151]	141	END SUBROUTINE check_conserve
[326]	142
	143	!------------------------ Detailed check (Lauritzen et al., 2014) -----------------------------
	144
[365]	145	SUBROUTINE check_conserve_detailed(it,tag, f_ps,f_dps,f_ue,f_theta_rhodz,f_phis)
[326]	146	USE pression_mod
	147	USE vorticity_mod
	148	USE caldyn_gcm_mod
	149	USE exner_mod
[365]	150	! USE mpipara, ONLY : is_mpi_root, comm_icosa
	151	USE omp_para, ONLY : is_master
	152	INTEGER, INTENT(IN) :: tag
[326]	153	TYPE(t_field),POINTER :: f_ps(:)
	154	TYPE(t_field),POINTER :: f_dps(:)
	155	TYPE(t_field),POINTER :: f_ue(:)
	156	TYPE(t_field),POINTER :: f_theta_rhodz(:)
	157	TYPE(t_field),POINTER :: f_phis(:)
	158	INTEGER, INTENT(IN) ::it
	159
	160	REAL(rstd),POINTER :: p(:,:),rhodz(:,:)
	161	INTEGER::ind,ierr
[365]	162	REAL(rstd) :: mtot, ztot, rmsdpdt, etot,stot,rmsv, ang_mass, ang_vel
[326]	163
	164	IF(check_type == check_detailed) THEN
	165
	166	CALL transfert_request(f_ue,req_e1_vect)
	167	CALL pression(f_ps,f_p)
	168	DO ind=1,ndomain
	169	IF (.NOT. assigned_domain(ind)) CYCLE
	170	CALL swap_dimensions(ind)
	171	CALL swap_geometry(ind)
	172	p=f_p(ind)
	173	rhodz=f_rhodz(ind)
	174	CALL compute_rhodz(p,rhodz)
	175	END DO
	176	CALL exner(f_ps,f_p,f_pks,f_pk)
[365]	177	CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vel, rmsv)
[326]	178
[365]	179	IF (is_master) THEN
	180	! XXX=mass,vel
	181	! AAM_XXX = value of total AAM_XXX at last check
	182	! ang_XXX = value of total AAM_XXX now
	183	! AAM_XXX_source(tag) = source of AAA_XXX due to tag=dissip,dyn,phys
	184	!WRITE(*,'(A,3E12.2)') 'AAM_mass', tag, ang_mass, AAM_mass, ang_mass-AAM_mass ! FIXME
	185	!WRITE(*,'(A,3E12.2)') 'AAM_vel ', tag, ang_vel, AAM_vel, ang_vel-AAM_vel ! FIXME
	186	AAM_mass_source(tag) = AAM_mass_source(tag) + ang_mass - AAM_mass
	187	AAM_vel_source(tag) = AAM_vel_source(tag) + ang_vel - AAM_vel
	188	AAM_mass = ang_mass
	189	AAM_vel = ang_vel
[326]	190	END IF
	191	END IF
	192	END SUBROUTINE check_conserve_detailed
	193
[149]	194	!---------------------------------------------------------------------
	195
[326]	196	SUBROUTINE global_sum(loc_sum, glob_sum)
	197	USE mpi_mod
	198	USE mpipara
	199	USE transfert_omp_mod
	200	REAL(rstd), INTENT(IN) :: loc_sum
	201	REAL(rstd), INTENT(OUT) :: glob_sum
	202	REAL(rstd) :: sum
	203
	204	CALL reduce_sum_omp(loc_sum, sum)
	205	IF (using_mpi) THEN
	206	!$OMP MASTER
	207	CALL MPI_REDUCE(sum, glob_sum, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr)
	208	!$OMP END MASTER
	209	ELSE
	210	glob_sum=loc_sum
[365]	211	ENDIF
[326]	212	END SUBROUTINE global_sum
	213
[160]	214	SUBROUTINE check_mass_conserve(f_ps,f_dps,mtot,rmsdpdt)
	215	USE mpi_mod
	216	USE mpipara
[186]	217	USE transfert_omp_mod
[295]	218	USE omp_para
[149]	219	IMPLICIT NONE
	220	TYPE(t_field),POINTER :: f_ps(:)
	221	TYPE(t_field),POINTER :: f_dps(:)
	222	REAL(rstd),POINTER :: ps(:),dps(:)
[160]	223	REAL(rstd), INTENT(OUT) :: mtot, rmsdpdt
	224
[149]	225	INTEGER :: ind,i,j,ij
[160]	226	REAL :: mloc, rmsloc
[186]	227	REAL :: mloc_mpi, rmsloc_mpi
[149]	228
[160]	229	mloc=0.0; rmsloc=0.0
[151]	230	DO ind=1,ndomain
[295]	231	IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
[149]	232	CALL swap_dimensions(ind)
	233	CALL swap_geometry(ind)
	234	dps=f_dps(ind)
	235	ps=f_ps(ind)
	236	DO j=jj_begin,jj_end
	237	DO i=ii_begin,ii_end
	238	ij=(j-1)*iim+i
	239	IF (domain(ind)%own(i,j)) THEN
[160]	240	mloc=mloc+ps(ij)*Ai(ij)
	241	rmsloc=rmsloc+dps(ij)*dps(ij)
[151]	242	ENDIF
[149]	243	ENDDO
	244	ENDDO
	245	ENDDO
[172]	246
[326]	247	CALL global_sum(mloc, mtot)
	248	CALL global_sum(rmsloc, rmsdpdt)
[149]	249	END SUBROUTINE check_mass_conserve
	250
	251	!---------------------------------------------------------------------
	252
[326]	253	SUBROUTINE check_energy(f_ue,f_theta_rhodz,f_phis, etot, &
	254	stot, AAM_mass_tot, AAM_vel_tot, rmsvtot)
[151]	255	USE pression_mod
	256	USE vorticity_mod
[149]	257	TYPE(t_field), POINTER :: f_ue(:)
	258	TYPE(t_field), POINTER :: f_theta_rhodz(:)
	259	TYPE(t_field), POINTER :: f_phis(:)
[326]	260	REAL(rstd),INTENT(OUT) :: etot, stot, AAM_mass_tot, AAM_vel_tot, rmsvtot
[149]	261	REAL(rstd), POINTER :: ue(:,:)
	262	REAL(rstd), POINTER :: p(:,:)
[387]	263	REAL(rstd), POINTER :: theta_rhodz(:,:,:)
[149]	264	REAL(rstd), POINTER :: pk(:,:)
	265	REAL(rstd), POINTER :: phis(:)
	266	REAL(rstd), POINTER :: rhodz(:,:)
[326]	267	REAL(rstd) :: e, s, AAM_mass, AAM_vel, rmsv
[149]	268	INTEGER :: ind
	269
[326]	270	e = 0 ; s = 0 ; AAM_mass = 0 ; AAM_vel=0 ; rmsv = 0
[149]	271
[151]	272	DO ind=1,ndomain
[186]	273	IF (.NOT. assigned_domain(ind)) CYCLE
[151]	274	CALL swap_dimensions(ind)
	275	CALL swap_geometry(ind)
	276	ue=f_ue(ind)
	277	p=f_p(ind)
	278	rhodz=f_rhodz(ind)
	279	theta_rhodz=f_theta_rhodz(ind)
	280	pk=f_pk(ind)
	281	phis=f_phis(ind)
[387]	282	CALL compute_energy(ind,ue,p,rhodz,theta_rhodz(:,:,1),pk,phis, &
[326]	283	e, s, AAM_mass, AAM_vel, rmsv)
[151]	284	END DO
[295]	285
[326]	286	CALL global_sum(e, etot)
	287	CALL global_sum(s, stot)
	288	CALL global_sum(rmsv, rmsvtot)
	289	CALL global_sum(AAM_mass, AAM_mass_tot)
	290	CALL global_sum(AAM_vel, AAM_vel_tot)
	291	END SUBROUTINE check_energy
[295]	292
[326]	293	SUBROUTINE compute_energy(ind,u,p,rhodz,theta_rhodz,pk,phis, e, s, AAM_mass, AAM_vel, rmsv)
	294	USE disvert_mod
	295	USE wind_mod
	296	USE omp_para
[149]	297	INTEGER,INTENT(IN)::ind
	298	REAL(rstd),INTENT(IN) :: u(3iimjjm,llm)
	299	REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
	300	REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm)
	301	REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm)
	302	REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm)
	303	REAL(rstd),INTENT(IN) :: phis(iim*jjm)
[326]	304	REAL(rstd),INTENT(INOUT) :: e, s, AAM_mass, AAM_vel, rmsv
[295]	305
[584]	306	REAL(rstd) :: ucenter(iim*jjm,llm,3)
[149]	307	REAL(rstd) :: ulon(iim*jjm,llm)
	308	REAL(rstd) :: ulat(iim*jjm,llm)
[326]	309
	310	REAL(rstd) :: mass_ij, theta_ij, KE_ij, rad,radomeg,lat
	311	INTEGER :: i,j,ij,l
	312
[151]	313	CALL compute_wind_centered(u,ucenter)
	314	CALL compute_wind_centered_lonlat_compound(ucenter,ulon,ulat)
	315	rad = radius
	316	radomeg = rad*omega
[326]	317
[295]	318	DO l = ll_begin,ll_end
[326]	319	DO j=jj_begin,jj_end
	320	DO i=ii_begin,ii_end
	321	ij=(j-1)*iim+i
	322	IF (domain(ind)%own(i,j)) THEN
	323	lat = lat_i(ij)
	324	mass_ij = rhodz(ij,l)*Ai(ij)
	325	theta_ij = theta_rhodz(ij,l)/rhodz(ij,l)
	326	KE_ij = 1/(4Ai(ij))(le(ij+u_right)de(ij+u_right)u(ij+u_right,l)**2 + &
	327	le(ij+u_rup)de(ij+u_rup)u(ij+u_rup,l)**2 + &
	328	le(ij+u_lup)de(ij+u_lup)u(ij+u_lup,l)**2 + &
	329	le(ij+u_left)de(ij+u_left)u(ij+u_left,l)**2 + &
	330	le(ij+u_ldown)de(ij+u_ldown)u(ij+u_ldown,l)**2 + &
	331	le(ij+u_rdown)de(ij+u_rdown)u(ij+u_rdown,l)**2 )
	332	rmsv = rmsv + 2mass_ijKE_ij
	333	s = s + Ai(ij)*theta_rhodz(ij,l)
	334	e = e + mass_ij(phis(ij)+theta_ijpk(ij,l)+KE_ij)
	335	AAM_mass = AAM_mass + radcos(lat)mass_ijradomegcos(lat)
	336	AAM_vel = AAM_vel + radcos(lat)mass_ij*ulon(ij,l)
	337	END IF
	338	END DO
	339	END DO
	340	END DO
	341
[365]	342	END SUBROUTINE compute_energy
[149]	343
	344	!---------------------------------------------------------------------
[365]	345
[295]	346	SUBROUTINE check_PV(ztot)
[365]	347	REAL(rstd),INTENT(OUT) :: ztot
	348	REAL(rstd), POINTER :: vort(:,:)
	349	REAL(rstd), POINTER :: rhodz(:,:)
	350	INTEGER :: ind
	351	REAL(rstd) :: z,z_mpi
	352
	353	z=0
	354	DO ind=1,ndomain
[186]	355	IF (.NOT. assigned_domain(ind)) CYCLE
[151]	356	CALL swap_dimensions(ind)
	357	CALL swap_geometry(ind)
	358	vort=f_vort(ind)
	359	rhodz=f_rhodz(ind)
[295]	360	CALL compute_PV(vort,rhodz,z)
[365]	361	ENDDO
	362	CALL global_sum(z, ztot)
[151]	363	END SUBROUTINE check_PV
[149]	364
[295]	365	SUBROUTINE compute_PV(vort,rhodz,z)
[151]	366	USE disvert_mod
[295]	367	USE omp_para
[149]	368	REAL(rstd),INTENT(IN) :: vort(2iimjjm,llm)
	369	REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm)
[295]	370	REAL(rstd),INTENT(INOUT) :: z
[149]	371	REAL(rstd)::qv1,qv2
	372	REAL(rstd)::hv1,hv2
	373	INTEGER :: i,j,ij,l,ij2
	374
[151]	375	hv1 = 0.0 ; hv2 = 0.0
[149]	376
[295]	377	DO l = ll_begin,ll_end
[149]	378	DO j=jj_begin+1,jj_end-1
[151]	379	DO i=ii_begin+1,ii_end-1
	380	ij=(j-1)*iim+i
	381	hv1 = Riv2(ij,vup) * rhodz(ij,l) &
	382	+ Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l) &
	383	+ Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
	384
	385	qv1 = ( vort(ij+z_up,l)+fv(ij+z_up) )/hv1
	386
	387	hv2 = Riv2(ij,vdown) * rhodz(ij,l) &
	388	+ Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l) &
	389	+ Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
	390
	391	qv2 =( vort(ij+z_down,l)+fv(ij+z_down) )/hv2
[149]	392
[295]	393	z = z + (qv1qv1hv1 + qv2qv2hv2)
[151]	394
[149]	395	ENDDO
[151]	396	ENDDO
[295]	397
[151]	398	ENDDO
	399
	400	END SUBROUTINE compute_PV
[149]	401	!---------------------------------------------------------------------
	402
[151]	403	SUBROUTINE compute_rhodz(p,rhodz)
[149]	404	REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
	405	REAL(rstd),INTENT(OUT):: rhodz(iim*jjm,llm)
	406	REAL(rstd) ::mass(iim*jjm,llm)
	407	INTEGER :: i,j,ij,l
	408
	409	DO l = 1, llm
	410	DO j=jj_begin,jj_end
	411	DO i=ii_begin,ii_end
	412	ij=(j-1)*iim+i
	413	mass(ij,l) = ( p(ij,l) - p(ij,l+1) ) * Ai(ij)/g
	414	rhodz(ij,l) = mass(ij,l) / Ai(ij)
	415	ENDDO
	416	ENDDO
	417	ENDDO
	418
[151]	419	END SUBROUTINE compute_rhodz
[149]	420	!====================================================================
	421
	422
	423	END MODULE check_conserve_mod

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

Download in other formats: