Context Navigation

source: codes/icosagcm/devel/src/dissip/dissip_gcm.f90 @ 614

Last change on this file since 614 was 550, checked in by dubos, 7 years ago
devel : backported commits 541 and 545 from trunk
File size: 26.0 KB

Rev	Line
[15]	1	MODULE dissip_gcm_mod
[19]	2	USE icosa
[15]	3
[26]	4	PRIVATE
	5
	6	TYPE(t_field),POINTER,SAVE :: f_due_diss1(:)
	7	TYPE(t_field),POINTER,SAVE :: f_due_diss2(:)
	8
	9	TYPE(t_field),POINTER,SAVE :: f_dtheta_diss(:)
	10	TYPE(t_field),POINTER,SAVE :: f_dtheta_rhodz_diss(:)
[186]	11	TYPE(t_message),SAVE :: req_due, req_dtheta
[15]	12
[26]	13	INTEGER,SAVE :: nitergdiv=1
[186]	14	!$OMP THREADPRIVATE(nitergdiv)
[26]	15	INTEGER,SAVE :: nitergrot=1
[186]	16	!$OMP THREADPRIVATE(nitergrot)
[26]	17	INTEGER,SAVE :: niterdivgrad=1
[186]	18	!$OMP THREADPRIVATE(niterdivgrad)
[15]	19
	20	REAL,ALLOCATABLE,SAVE :: tau_graddiv(:)
[186]	21	!$OMP THREADPRIVATE(tau_graddiv)
[15]	22	REAL,ALLOCATABLE,SAVE :: tau_gradrot(:)
[186]	23	!$OMP THREADPRIVATE(tau_gradrot)
[15]	24	REAL,ALLOCATABLE,SAVE :: tau_divgrad(:)
[186]	25	!$OMP THREADPRIVATE(tau_divgrad)
[15]	26
	27	REAL,SAVE :: cgraddiv
[186]	28	!$OMP THREADPRIVATE(cgraddiv)
[15]	29	REAL,SAVE :: cgradrot
[186]	30	!$OMP THREADPRIVATE(cgradrot)
[15]	31	REAL,SAVE :: cdivgrad
[186]	32	!$OMP THREADPRIVATE(cdivgrad)
[109]	33
	34	INTEGER, SAVE :: rayleigh_friction_type, rayleigh_shear
[186]	35	!$OMP THREADPRIVATE(rayleigh_friction_type)
[529]	36	!$OMP THREADPRIVATE(rayleigh_shear)
[186]	37	REAL, SAVE :: rayleigh_tau
[529]	38	!$OMP THREADPRIVATE(rayleigh_tau)
[15]	39
[26]	40	REAL,SAVE :: dtdissip
[186]	41	!$OMP THREADPRIVATE(dtdissip)
[151]	42
[26]	43	PUBLIC init_dissip, dissip
[15]	44	CONTAINS
	45
	46	SUBROUTINE allocate_dissip
[19]	47	USE icosa
[15]	48	IMPLICIT NONE
[26]	49	CALL allocate_field(f_due_diss1,field_u,type_real,llm)
	50	CALL allocate_field(f_due_diss2,field_u,type_real,llm)
	51	CALL allocate_field(f_dtheta_diss,field_t,type_real,llm)
	52	CALL allocate_field(f_dtheta_rhodz_diss,field_t,type_real,llm)
[15]	53	ALLOCATE(tau_graddiv(llm))
	54	ALLOCATE(tau_gradrot(llm))
	55	ALLOCATE(tau_divgrad(llm))
	56	END SUBROUTINE allocate_dissip
	57
[98]	58	SUBROUTINE init_dissip
[19]	59	USE icosa
[15]	60	USE disvert_mod
[26]	61	USE mpi_mod
	62	USE mpipara
[151]	63	USE transfert_mod
[148]	64	USE time_mod
[186]	65	USE transfert_omp_mod
[295]	66	USE omp_para
[15]	67	IMPLICIT NONE
	68
[186]	69	TYPE(t_field),POINTER,SAVE :: f_u(:)
	70	TYPE(t_field),POINTER,SAVE :: f_du(:)
	71	REAL(rstd),POINTER :: u(:)
	72	REAL(rstd),POINTER :: du(:)
	73	TYPE(t_field),POINTER,SAVE :: f_theta(:)
	74	TYPE(t_field),POINTER ,SAVE :: f_dtheta(:)
[15]	75	REAL(rstd),POINTER :: theta(:)
	76	REAL(rstd),POINTER :: dtheta(:)
[59]	77	REAL(rstd) :: dumax,dumax1
	78	REAL(rstd) :: dthetamax,dthetamax1
[15]	79	REAL(rstd) :: r
	80	REAL(rstd) :: tau
	81	REAL(rstd) :: zz, zvert, fact
	82	INTEGER :: l
[109]	83	CHARACTER(len=255) :: rayleigh_friction_key
[148]	84	REAL(rstd) :: mintau
[174]	85	INTEGER :: seed_size
	86	INTEGER,ALLOCATABLE :: seed(:)
[15]	87
	88
[550]	89	INTEGER :: i,j,ij,ind,it,iter,M
[15]	90
[109]	91	rayleigh_friction_key='none'
	92	CALL getin("rayleigh_friction_type",rayleigh_friction_key)
	93	SELECT CASE(TRIM(rayleigh_friction_key))
	94	CASE('none')
	95	rayleigh_friction_type=0
[295]	96	IF (is_master) PRINT *, 'No Rayleigh friction'
[109]	97	CASE('dcmip2_schaer_noshear')
	98	rayleigh_friction_type=1
	99	rayleigh_shear=0
[295]	100	IF (is_master) PRINT *, 'Rayleigh friction : Schaer-like mountain without shear DCMIP2.1'
[109]	101	CASE('dcmip2_schaer_shear')
	102	rayleigh_shear=1
	103	rayleigh_friction_type=2
[295]	104	IF (is_master) PRINT *, 'Rayleigh friction : Schaer-like mountain with shear DCMIP2.2'
[109]	105	CASE DEFAULT
[358]	106	IF (is_master) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), &
	107	' in dissip_gcm.f90/init_dissip'
[109]	108	STOP
	109	END SELECT
	110
	111	IF(rayleigh_friction_type>0) THEN
	112	rayleigh_tau=0.
	113	CALL getin("rayleigh_friction_tau",rayleigh_tau)
	114	rayleigh_tau = rayleigh_tau / scale_factor
	115	IF(rayleigh_tau<=0) THEN
[295]	116	IF (is_master) PRINT *, 'Forbidden : negative value for rayleigh_friction_tau =',rayleigh_tau
[109]	117	STOP
	118	END IF
	119	END IF
	120
[15]	121	CALL allocate_dissip
	122	CALL allocate_field(f_u,field_u,type_real)
	123	CALL allocate_field(f_du,field_u,type_real)
	124	CALL allocate_field(f_theta,field_t,type_real)
	125	CALL allocate_field(f_dtheta,field_t,type_real)
[151]	126
	127	CALL init_message(f_due_diss1,req_e1_vect,req_due)
	128	CALL init_message(f_dtheta_diss,req_i1,req_dtheta)
[15]	129
	130	tau_graddiv(:)=5000
	131	CALL getin("tau_graddiv",tau)
[32]	132	tau_graddiv(:)=tau/scale_factor
[26]	133
	134	CALL getin("nitergdiv",nitergdiv)
[15]	135
	136	tau_gradrot(:)=5000
[186]	137	CALL getin("tau_gradrot",tau)
[32]	138	tau_gradrot(:)=tau/scale_factor
[26]	139
	140	CALL getin("nitergrot",nitergrot)
[15]	141
	142
	143	tau_divgrad(:)=5000
	144	CALL getin("tau_divgrad",tau)
[32]	145	tau_divgrad(:)=tau/scale_factor
[26]	146
	147	CALL getin("niterdivgrad",niterdivgrad)
[15]	148
	149
	150	cgraddiv=-1
	151	cdivgrad=-1
	152	cgradrot=-1
[295]	153
	154	!$OMP BARRIER
	155	!$OMP MASTER
[15]	156	DO ind=1,ndomain
	157	CALL swap_dimensions(ind)
	158	CALL swap_geometry(ind)
	159	u=f_u(ind)
	160
[550]	161	! set random seed to get reproductibility when using a different number of process
	162	CALL RANDOM_SEED(size=M)
	163	CALL RANDOM_SEED(put=(/(i,i=1,M)/))
	164
[15]	165	DO j=jj_begin,jj_end
	166	DO i=ii_begin,ii_end
[174]	167	ij=(j-1)*iim+i
[15]	168	CALL RANDOM_NUMBER(r)
[174]	169	u(ij+u_right)=r-0.5
[15]	170	CALL RANDOM_NUMBER(r)
[174]	171	u(ij+u_lup)=r-0.5
[15]	172	CALL RANDOM_NUMBER(r)
[174]	173	u(ij+u_ldown)=r-0.5
	174	ENDDO
	175	ENDDO
[15]	176	ENDDO
[295]	177	!$OMP END MASTER
	178	!$OMP BARRIER
[15]	179
[29]	180	DO it=1,20
[26]	181
[15]	182	dumax=0
[26]	183	DO iter=1,nitergdiv
[146]	184	CALL transfert_request(f_u,req_e1_vect)
[26]	185	DO ind=1,ndomain
[295]	186	IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
[26]	187	CALL swap_dimensions(ind)
	188	CALL swap_geometry(ind)
	189	u=f_u(ind)
	190	du=f_du(ind)
[151]	191	CALL compute_gradiv(u,du,1,1)
[26]	192	u=du
	193	ENDDO
[15]	194	ENDDO
[26]	195
[146]	196	CALL transfert_request(f_du,req_e1_vect)
[15]	197
	198	DO ind=1,ndomain
[295]	199	IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
[15]	200	CALL swap_dimensions(ind)
	201	CALL swap_geometry(ind)
	202	u=f_u(ind)
	203	du=f_du(ind)
[26]	204
[15]	205	DO j=jj_begin,jj_end
	206	DO i=ii_begin,ii_end
[174]	207	ij=(j-1)*iim+i
	208	if (le(ij+u_right)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_right)))
	209	if (le(ij+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_lup)))
	210	if (le(ij+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_ldown)))
[15]	211	ENDDO
	212	ENDDO
	213	ENDDO
	214
[59]	215	IF (using_mpi) THEN
[295]	216	CALL reduce_max_omp(dumax,dumax1)
[186]	217	!$OMP MASTER
	218	CALL MPI_ALLREDUCE(dumax1,dumax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
	219	!$OMP END MASTER
	220	CALL bcast_omp(dumax)
	221	ELSE
[295]	222	CALL allreduce_max_omp(dumax,dumax1)
[59]	223	dumax=dumax1
[186]	224	ENDIF
[59]	225
[15]	226	DO ind=1,ndomain
[295]	227	IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
[15]	228	CALL swap_dimensions(ind)
	229	CALL swap_geometry(ind)
	230	u=f_u(ind)
	231	du=f_du(ind)
	232	u=du/dumax
	233	ENDDO
[295]	234	IF (is_master) PRINT *,"gradiv : it :",it ,": dumax",dumax
[15]	235
	236	ENDDO
[295]	237	IF (is_master) PRINT *,"gradiv : dumax",dumax
	238	IF (is_master) PRINT , 'mean T-cell edge size (km)', 1.45radius/iim_glo/1000., &
[131]	239	'effective T-cell half-edge size (km)', dumax**(-.5/nitergdiv)/1000
[387]	240	IF (is_master) PRINT *, 'Max. time step assuming c=340 m/s and Courant number=3 (ARK2.3) :', &
	241	3./340.dumax*(-.5/nitergdiv)
[129]	242
[15]	243	cgraddiv=dumax**(-1./nitergdiv)
[295]	244	IF (is_master) PRINT *,"cgraddiv : ",cgraddiv
[15]	245
[295]	246	!$OMP BARRIER
	247	!$OMP MASTER
[15]	248	DO ind=1,ndomain
	249	CALL swap_dimensions(ind)
	250	CALL swap_geometry(ind)
	251	u=f_u(ind)
	252
[550]	253	! set random seed to get reproductibility when using a different number of process
	254	CALL RANDOM_SEED(size=M)
	255	CALL RANDOM_SEED(put=(/(i,i=1,M)/))
	256
	257	DO j=jj_begin,jj_end
[15]	258	DO i=ii_begin,ii_end
[174]	259	ij=(j-1)*iim+i
[15]	260	CALL RANDOM_NUMBER(r)
[174]	261	u(ij+u_right)=r-0.5
[15]	262	CALL RANDOM_NUMBER(r)
[174]	263	u(ij+u_lup)=r-0.5
[15]	264	CALL RANDOM_NUMBER(r)
[174]	265	u(ij+u_ldown)=r-0.5
	266	ENDDO
	267	ENDDO
[15]	268	ENDDO
[295]	269	!$OMP END MASTER
	270	!$OMP BARRIER
[15]	271
	272
[29]	273	DO it=1,20
[15]	274
	275	dumax=0
[26]	276	DO iter=1,nitergrot
[146]	277	CALL transfert_request(f_u,req_e1_vect)
[26]	278	DO ind=1,ndomain
[295]	279	IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
[26]	280	CALL swap_dimensions(ind)
	281	CALL swap_geometry(ind)
	282	u=f_u(ind)
	283	du=f_du(ind)
[151]	284	CALL compute_gradrot(u,du,1,1)
[26]	285	u=du
	286	ENDDO
[15]	287	ENDDO
	288
[146]	289	CALL transfert_request(f_du,req_e1_vect)
[15]	290
	291	DO ind=1,ndomain
[295]	292	IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
[15]	293	CALL swap_dimensions(ind)
	294	CALL swap_geometry(ind)
	295	u=f_u(ind)
	296	du=f_du(ind)
	297
	298	DO j=jj_begin,jj_end
	299	DO i=ii_begin,ii_end
[174]	300	ij=(j-1)*iim+i
	301	if (le(ij+u_right)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_right)))
	302	if (le(ij+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_lup)))
	303	if (le(ij+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_ldown)))
[15]	304	ENDDO
	305	ENDDO
	306	ENDDO
[26]	307
[186]	308	IF (using_mpi) THEN
[295]	309	CALL reduce_max_omp(dumax,dumax1)
[186]	310	!$OMP MASTER
	311	CALL MPI_ALLREDUCE(dumax1,dumax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
	312	!$OMP END MASTER
	313	CALL bcast_omp(dumax)
	314	ELSE
[295]	315	CALL allreduce_max_omp(dumax,dumax1)
[59]	316	dumax=dumax1
[186]	317	ENDIF
	318
[15]	319
	320	DO ind=1,ndomain
[295]	321	IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
[15]	322	CALL swap_dimensions(ind)
	323	CALL swap_geometry(ind)
	324	u=f_u(ind)
	325	du=f_du(ind)
	326	u=du/dumax
	327	ENDDO
	328
[295]	329	IF (is_master) PRINT *,"gradrot : it :",it ,": dumax",dumax
[15]	330
	331	ENDDO
[295]	332	IF (is_master) PRINT *,"gradrot : dumax",dumax
[15]	333
[26]	334	cgradrot=dumax**(-1./nitergrot)
[295]	335	IF (is_master) PRINT *,"cgradrot : ",cgradrot
[15]	336
	337
	338
[295]	339	!$OMP BARRIER
	340	!$OMP MASTER
[15]	341	DO ind=1,ndomain
	342	CALL swap_dimensions(ind)
	343	CALL swap_geometry(ind)
	344	theta=f_theta(ind)
[550]	345
	346	! set random seed to get reproductibility when using a different number of process
	347	CALL RANDOM_SEED(size=M)
	348	CALL RANDOM_SEED(put=(/(i,i=1,M)/))
	349
[174]	350	DO j=jj_begin,jj_end
	351	DO i=ii_begin,ii_end
	352	ij=(j-1)*iim+i
	353	CALL RANDOM_NUMBER(r)
	354	theta(ij)=r-0.5
[15]	355	ENDDO
[174]	356	ENDDO
[15]	357	ENDDO
[295]	358	!$OMP END MASTER
	359	!$OMP BARRIER
[15]	360
[29]	361	DO it=1,20
[15]	362
	363	dthetamax=0
[26]	364	DO iter=1,niterdivgrad
	365	CALL transfert_request(f_theta,req_i1)
	366	DO ind=1,ndomain
[295]	367	IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
[26]	368	CALL swap_dimensions(ind)
	369	CALL swap_geometry(ind)
	370	theta=f_theta(ind)
	371	dtheta=f_dtheta(ind)
[151]	372	CALL compute_divgrad(theta,dtheta,1,1)
[26]	373	theta=dtheta
	374	ENDDO
[15]	375	ENDDO
	376
	377	CALL transfert_request(f_dtheta,req_i1)
	378
	379	DO ind=1,ndomain
[295]	380	IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
[15]	381	CALL swap_dimensions(ind)
	382	CALL swap_geometry(ind)
	383	theta=f_theta(ind)
	384	dtheta=f_dtheta(ind)
	385
	386	DO j=jj_begin,jj_end
	387	DO i=ii_begin,ii_end
[174]	388	ij=(j-1)*iim+i
	389	dthetamax=MAX(dthetamax,ABS(dtheta(ij)))
[15]	390	ENDDO
	391	ENDDO
	392	ENDDO
[186]	393
	394	IF (using_mpi) THEN
[295]	395	CALL reduce_max_omp(dthetamax ,dthetamax1)
[186]	396	!$OMP MASTER
	397	CALL MPI_ALLREDUCE(dthetamax1,dthetamax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
	398	!$OMP END MASTER
	399	CALL bcast_omp(dthetamax)
	400	ELSE
[295]	401	CALL allreduce_max_omp(dthetamax,dthetamax1)
[186]	402	dumax=dumax1
	403	ENDIF
[174]	404
[295]	405	IF (is_master) PRINT *,"divgrad : it :",it ,": dthetamax",dthetamax
[15]	406
	407	DO ind=1,ndomain
[295]	408	IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
[15]	409	CALL swap_dimensions(ind)
	410	CALL swap_geometry(ind)
	411	theta=f_theta(ind)
	412	dtheta=f_dtheta(ind)
	413	theta=dtheta/dthetamax
	414	ENDDO
	415	ENDDO
	416
[295]	417	IF (is_master) PRINT *,"divgrad : divgrad",dthetamax
[15]	418
[26]	419	cdivgrad=dthetamax**(-1./niterdivgrad)
[295]	420	IF (is_master) PRINT *,"cdivgrad : ",cdivgrad
[15]	421
	422
	423	fact=2
	424	DO l=1,llm
[167]	425	IF(ap_bp_present) THEN ! height-dependent dissipation
	426	zz= 1. - preff/presnivs(l)
	427	ELSE
	428	zz = 0.
	429	END IF
[15]	430	zvert=fact-(fact-1)/(1+zz*zz)
	431	tau_graddiv(l) = tau_graddiv(l)/zvert
	432	tau_gradrot(l) = tau_gradrot(l)/zvert
	433	tau_divgrad(l) = tau_divgrad(l)/zvert
	434	ENDDO
[148]	435
	436	mintau=tau_graddiv(1)
	437	DO l=1,llm
	438	mintau=MIN(mintau,tau_graddiv(l))
	439	mintau=MIN(mintau,tau_gradrot(l))
	440	mintau=MIN(mintau,tau_divgrad(l))
	441	ENDDO
[529]	442
	443	IF(rayleigh_friction_type>0) mintau=MIN(mintau,rayleigh_tau)
	444
[250]	445	IF(mintau>0) THEN
	446	itau_dissip=INT(mintau/dt)
	447	dtdissip=itau_dissip*dt
	448	ELSE
[295]	449	IF (is_master) PRINT *,"No dissipation time set, setting itau_dissip to 1000000000"
[250]	450	itau_dissip=100000000
	451	END IF
[148]	452	itau_dissip=MAX(1,itau_dissip)
[529]	453	IF (is_master) PRINT *,"rayleigh_tau",rayleigh_tau, "mintau ",mintau, &
	454	"itau_dissip",itau_dissip," dtdissip ",dtdissip
[15]	455
	456	END SUBROUTINE init_dissip
	457
	458
[529]	459	SUBROUTINE dissip(f_ps,f_mass,f_phis,f_geopot,f_theta_rhodz,f_ue, f_dtheta_rhodz,f_due)
[19]	460	USE icosa
[15]	461	USE theta2theta_rhodz_mod
[109]	462	USE pression_mod
	463	USE exner_mod
	464	USE geopotential_mod
[145]	465	USE trace
[148]	466	USE time_mod
[151]	467	USE omp_para
[15]	468	IMPLICIT NONE
[529]	469	TYPE(t_field),POINTER :: f_ps(:), f_mass(:), f_phis(:), f_geopot(:)
	470	TYPE(t_field),POINTER :: f_theta_rhodz(:), f_dtheta_rhodz(:)
	471	TYPE(t_field),POINTER :: f_ue(:), f_due(:)
[26]	472
[15]	473	REAL(rstd),POINTER :: due(:,:)
[109]	474	REAL(rstd),POINTER :: phi(:,:), ue(:,:)
[26]	475	REAL(rstd),POINTER :: due_diss1(:,:)
	476	REAL(rstd),POINTER :: due_diss2(:,:)
[387]	477	REAL(rstd),POINTER :: dtheta_rhodz(:,:,:)
[26]	478	REAL(rstd),POINTER :: dtheta_rhodz_diss(:,:)
[15]	479
[109]	480	INTEGER :: ind, shear
[174]	481	INTEGER :: l,ij
[151]	482
	483	!$OMP BARRIER
[15]	484
[145]	485	CALL trace_start("dissip")
[26]	486	CALL gradiv(f_ue,f_due_diss1)
	487	CALL gradrot(f_ue,f_due_diss2)
[167]	488
	489	CALL divgrad_theta_rhodz(f_mass,f_theta_rhodz,f_dtheta_rhodz_diss)
[26]	490
[15]	491	DO ind=1,ndomain
[186]	492	IF (.NOT. assigned_domain(ind)) CYCLE
[15]	493	CALL swap_dimensions(ind)
	494	CALL swap_geometry(ind)
	495	due=f_due(ind)
[26]	496	due_diss1=f_due_diss1(ind)
	497	due_diss2=f_due_diss2(ind)
	498	dtheta_rhodz=f_dtheta_rhodz(ind)
	499	dtheta_rhodz_diss=f_dtheta_rhodz_diss(ind)
[151]	500
	501	DO l=ll_begin,ll_end
[487]	502	!DIR$ SIMD
[174]	503	DO ij=ij_begin,ij_end
[15]	504
[174]	505	due(ij+u_right,l) = -0.5( due_diss1(ij+u_right,l)/tau_graddiv(l) + due_diss2(ij+u_right,l)/tau_gradrot(l))itau_dissip
	506	due(ij+u_lup,l) = -0.5( due_diss1(ij+u_lup,l) /tau_graddiv(l) + due_diss2(ij+u_lup,l) /tau_gradrot(l))itau_dissip
	507	due(ij+u_ldown,l) = -0.5( due_diss1(ij+u_ldown,l)/tau_graddiv(l) + due_diss2(ij+u_ldown,l)/tau_gradrot(l))itau_dissip
[26]	508
[387]	509	dtheta_rhodz(ij,l,1) = -0.5dtheta_rhodz_diss(ij,l)/tau_divgrad(l)itau_dissip
[26]	510	ENDDO
	511	ENDDO
[109]	512
[151]	513	! dtheta_rhodz=0
	514	! due=0
	515
[529]	516	IF(rayleigh_friction_type>0) THEN
	517	phi=f_geopot(ind)
	518	ue=f_ue(ind)
	519	DO l=ll_begin,ll_end
	520	DO ij=ij_begin,ij_end
	521	CALL relax(t_right, u_right)
	522	CALL relax(t_lup, u_lup)
	523	CALL relax(t_ldown, u_ldown)
	524	ENDDO
	525	END DO
	526	END IF
[109]	527	END DO
	528
[145]	529	CALL trace_end("dissip")
	530
[151]	531	!$OMP BARRIER
	532
[109]	533	CONTAINS
	534	SUBROUTINE relax(shift_t, shift_u)
	535	USE dcmip_initial_conditions_test_1_2_3
	536	REAL(rstd) :: z, ulon,ulat, lon,lat, & ! input to test2_schaer_mountain
	537	p,hyam,hybm,w,t,phis,ps,rho,q, & ! unused input/output to test2_schaer_mountain
	538	fz, u3d(3), uref
	539	REAL(rstd), PARAMETER :: zh=2e4,ztop=3e4 ! DCMIP values
	540	LOGICAL :: hybrid_eta
	541	INTEGER :: shift_u, shift_t, zcoords, nn
[174]	542	z = (phi(ij,l)+phi(ij+shift_t,l))/(2.*g)
[109]	543	IF(z>zh) THEN ! relax only in the sponge layer z>zh
[151]	544
[174]	545	nn = ij+shift_u
[109]	546	zcoords = 1 ; hybrid_eta = .FALSE. ! use z instead of p or hyam/hybm
[286]	547	CALL test2_schaer_mountain(lon_e(nn),lat_e(nn),p,z,zcoords,hybrid_eta, &
[109]	548	hyam,hybm,shear,ulon,ulat,w,t,phis,ps,rho,q)
	549	! u3d = ulonelon_e(nn,:) + ulatelat_e(nn,:)
	550	u3d = ulon*elon_e(nn,:) ! ulat=0
	551	uref = sum(u3d*ep_e(nn,:))
	552
	553	fz = sin((pi/2)*(z-zh)/(ztop-zh))
	554	fz = fz*fz/rayleigh_tau
[529]	555	due(nn,l) = due(nn,l) - itau_dissipfz(ue(nn,l)-uref)
	556	! fz = 1./rayleigh_tau
[109]	557	! due(nn,l) = due(nn,l) - fz*ue(nn,l)
	558	END IF
	559	END SUBROUTINE relax
[15]	560
	561	END SUBROUTINE dissip
	562
[26]	563	SUBROUTINE gradiv(f_ue,f_due)
[19]	564	USE icosa
[145]	565	USE trace
[151]	566	USE omp_para
[15]	567	IMPLICIT NONE
[26]	568	TYPE(t_field),POINTER :: f_ue(:)
	569	TYPE(t_field),POINTER :: f_due(:)
	570	REAL(rstd),POINTER :: ue(:,:)
	571	REAL(rstd),POINTER :: due(:,:)
	572	INTEGER :: ind
[174]	573	INTEGER :: it,l,ij
[26]	574
[145]	575	CALL trace_start("gradiv")
	576
[26]	577	DO ind=1,ndomain
[186]	578	IF (.NOT. assigned_domain(ind)) CYCLE
[26]	579	CALL swap_dimensions(ind)
	580	CALL swap_geometry(ind)
	581	ue=f_ue(ind)
	582	due=f_due(ind)
[151]	583	DO l = ll_begin, ll_end
[487]	584	!DIR$ SIMD
[174]	585	DO ij=ij_begin,ij_end
[151]	586	due(ij+u_right,l)=ue(ij+u_right,l)
	587	due(ij+u_lup,l)=ue(ij+u_lup,l)
	588	due(ij+u_ldown,l)=ue(ij+u_ldown,l)
	589	ENDDO
	590	ENDDO
[26]	591	ENDDO
[15]	592
[26]	593	DO it=1,nitergdiv
	594
[151]	595	CALL send_message(f_due,req_due)
	596	CALL wait_message(req_due)
[26]	597
	598	DO ind=1,ndomain
[186]	599	IF (.NOT. assigned_domain(ind)) CYCLE
[26]	600	CALL swap_dimensions(ind)
	601	CALL swap_geometry(ind)
	602	due=f_due(ind)
[151]	603	CALL compute_gradiv(due,due,ll_begin,ll_end)
[26]	604	ENDDO
	605	ENDDO
[15]	606
[145]	607	CALL trace_end("gradiv")
	608
[26]	609	END SUBROUTINE gradiv
	610
	611
	612	SUBROUTINE gradrot(f_ue,f_due)
	613	USE icosa
[145]	614	USE trace
[151]	615	USE omp_para
[26]	616	IMPLICIT NONE
	617	TYPE(t_field),POINTER :: f_ue(:)
	618	TYPE(t_field),POINTER :: f_due(:)
	619	REAL(rstd),POINTER :: ue(:,:)
	620	REAL(rstd),POINTER :: due(:,:)
[15]	621	INTEGER :: ind
[174]	622	INTEGER :: it,l,ij
[26]	623
[145]	624	CALL trace_start("gradrot")
	625
[26]	626	DO ind=1,ndomain
[186]	627	IF (.NOT. assigned_domain(ind)) CYCLE
[26]	628	CALL swap_dimensions(ind)
	629	CALL swap_geometry(ind)
	630	ue=f_ue(ind)
	631	due=f_due(ind)
[151]	632	DO l = ll_begin, ll_end
[487]	633	!DIR$ SIMD
[174]	634	DO ij=ij_begin,ij_end
[151]	635	due(ij+u_right,l)=ue(ij+u_right,l)
	636	due(ij+u_lup,l)=ue(ij+u_lup,l)
	637	due(ij+u_ldown,l)=ue(ij+u_ldown,l)
	638	ENDDO
	639	ENDDO
[26]	640	ENDDO
[15]	641
[26]	642	DO it=1,nitergrot
	643
[151]	644	CALL send_message(f_due,req_due)
	645	CALL wait_message(req_due)
[26]	646
	647	DO ind=1,ndomain
[186]	648	IF (.NOT. assigned_domain(ind)) CYCLE
[26]	649	CALL swap_dimensions(ind)
	650	CALL swap_geometry(ind)
	651	due=f_due(ind)
[151]	652	CALL compute_gradrot(due,due,ll_begin,ll_end)
[15]	653	ENDDO
	654
[26]	655	ENDDO
[15]	656
[145]	657	CALL trace_end("gradrot")
	658
[26]	659	END SUBROUTINE gradrot
	660
	661	SUBROUTINE divgrad(f_theta,f_dtheta)
	662	USE icosa
[145]	663	USE trace
[151]	664	USE omp_para
[26]	665	IMPLICIT NONE
	666	TYPE(t_field),POINTER :: f_theta(:)
	667	TYPE(t_field),POINTER :: f_dtheta(:)
	668	REAL(rstd),POINTER :: theta(:,:)
	669	REAL(rstd),POINTER :: dtheta(:,:)
	670	INTEGER :: ind
	671	INTEGER :: it
[145]	672
	673	CALL trace_start("divgrad")
[26]	674
	675	DO ind=1,ndomain
[186]	676	IF (.NOT. assigned_domain(ind)) CYCLE
[26]	677	CALL swap_dimensions(ind)
	678	CALL swap_geometry(ind)
	679	theta=f_theta(ind)
	680	dtheta=f_dtheta(ind)
	681	dtheta=theta
	682	ENDDO
[15]	683
[26]	684	DO it=1,niterdivgrad
	685
	686	CALL transfert_request(f_dtheta,req_i1)
	687
	688	DO ind=1,ndomain
[186]	689	IF (.NOT. assigned_domain(ind)) CYCLE
[26]	690	CALL swap_dimensions(ind)
	691	CALL swap_geometry(ind)
	692	dtheta=f_dtheta(ind)
[151]	693	CALL compute_divgrad(dtheta,dtheta,ll_begin,ll_end)
[26]	694	ENDDO
[15]	695
[26]	696	ENDDO
	697
[145]	698	CALL trace_end("divgrad")
	699
[26]	700	END SUBROUTINE divgrad
	701
[167]	702	SUBROUTINE divgrad_theta_rhodz(f_mass,f_theta_rhodz,f_dtheta_rhodz)
[151]	703	USE icosa
	704	USE trace
	705	USE omp_para
	706	IMPLICIT NONE
[167]	707	TYPE(t_field),POINTER :: f_mass(:)
[151]	708	TYPE(t_field),POINTER :: f_theta_rhodz(:)
	709	TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
	710
[167]	711	REAL(rstd),POINTER :: mass(:,:)
[387]	712	REAL(rstd),POINTER :: theta_rhodz(:,:,:)
[151]	713	REAL(rstd),POINTER :: dtheta_rhodz(:,:)
[26]	714
[151]	715	INTEGER :: ind
[174]	716	INTEGER :: it,l,ij
[151]	717
	718	CALL trace_start("divgrad")
[295]	719
[151]	720	DO ind=1,ndomain
[186]	721	IF (.NOT. assigned_domain(ind)) CYCLE
[151]	722	CALL swap_dimensions(ind)
	723	CALL swap_geometry(ind)
[167]	724	mass=f_mass(ind)
[151]	725	theta_rhodz=f_theta_rhodz(ind)
	726	dtheta_rhodz=f_dtheta_rhodz(ind)
	727	DO l = ll_begin, ll_end
[487]	728	!DIR$ SIMD
[174]	729	DO ij=ij_begin,ij_end
[387]	730	dtheta_rhodz(ij,l) = theta_rhodz(ij,l,1) / mass(ij,l)
[151]	731	ENDDO
	732	ENDDO
	733	ENDDO
	734
	735	DO it=1,niterdivgrad
	736
	737	CALL send_message(f_dtheta_rhodz,req_dtheta)
	738	CALL wait_message(req_dtheta)
	739	DO ind=1,ndomain
[186]	740	IF (.NOT. assigned_domain(ind)) CYCLE
[151]	741	CALL swap_dimensions(ind)
	742	CALL swap_geometry(ind)
	743	dtheta_rhodz=f_dtheta_rhodz(ind)
	744	CALL compute_divgrad(dtheta_rhodz,dtheta_rhodz,ll_begin,ll_end)
	745	ENDDO
	746
	747	ENDDO
	748
	749	DO ind=1,ndomain
[186]	750	IF (.NOT. assigned_domain(ind)) CYCLE
[151]	751	CALL swap_dimensions(ind)
	752	CALL swap_geometry(ind)
	753	dtheta_rhodz=f_dtheta_rhodz(ind)
[175]	754	mass=f_mass(ind)
[151]	755
	756	DO l = ll_begin, ll_end
[487]	757	!DIR$ SIMD
[174]	758	DO ij=ij_begin,ij_end
[167]	759	dtheta_rhodz(ij,l) = dtheta_rhodz(ij,l) * mass(ij,l)
[151]	760	ENDDO
	761	ENDDO
	762	ENDDO
	763
	764
	765	CALL trace_end("divgrad")
	766
	767	END SUBROUTINE divgrad_theta_rhodz
[15]	768
[151]	769	SUBROUTINE compute_gradiv(ue,gradivu_e,llb,lle)
[19]	770	USE icosa
[15]	771	IMPLICIT NONE
[151]	772	INTEGER,INTENT(IN) :: llb
	773	INTEGER,INTENT(IN) :: lle
	774	REAL(rstd),INTENT(IN) :: ue(iim3jjm,llm)
	775	REAL(rstd),INTENT(OUT) :: gradivu_e(iim3jjm,llm)
	776	REAL(rstd) :: divu_i(iim*jjm,llb:lle)
[15]	777
[174]	778	INTEGER :: ij,l
[15]	779
[151]	780	DO l=llb,lle
[487]	781	!DIR$ SIMD
[174]	782	DO ij=ij_begin,ij_end
	783	divu_i(ij,l)=1./Ai(ij)(ne(ij,right)ue(ij+u_right,l)*le(ij+u_right) + &
	784	ne(ij,rup)ue(ij+u_rup,l)le(ij+u_rup) + &
	785	ne(ij,lup)ue(ij+u_lup,l)le(ij+u_lup) + &
	786	ne(ij,left)ue(ij+u_left,l)le(ij+u_left) + &
	787	ne(ij,ldown)ue(ij+u_ldown,l)le(ij+u_ldown) + &
	788	ne(ij,rdown)ue(ij+u_rdown,l)le(ij+u_rdown))
[15]	789	ENDDO
	790	ENDDO
	791
[151]	792	DO l=llb,lle
[487]	793	!DIR$ SIMD
[174]	794	DO ij=ij_begin,ij_end
[15]	795
[174]	796	gradivu_e(ij+u_right,l)=-1/de(ij+u_right)(ne(ij,right)divu_i(ij,l)+ ne(ij+t_right,left)*divu_i(ij+t_right,l) )
[15]	797
[174]	798	gradivu_e(ij+u_lup,l)=-1/de(ij+u_lup)(ne(ij,lup)divu_i(ij,l)+ ne(ij+t_lup,rdown)*divu_i(ij+t_lup,l))
[15]	799
[174]	800	gradivu_e(ij+u_ldown,l)=-1/de(ij+u_ldown)(ne(ij,ldown)divu_i(ij,l)+ne(ij+t_ldown,rup)*divu_i(ij+t_ldown,l) )
[15]	801
	802	ENDDO
	803	ENDDO
	804
[151]	805	DO l=llb,lle
[487]	806	!DIR$ SIMD
[174]	807	DO ij=ij_begin,ij_end
	808	gradivu_e(ij+u_right,l)=-gradivu_e(ij+u_right,l)*cgraddiv
	809	gradivu_e(ij+u_lup,l)=-gradivu_e(ij+u_lup,l)*cgraddiv
	810	gradivu_e(ij+u_ldown,l)=-gradivu_e(ij+u_ldown,l)*cgraddiv
[15]	811	ENDDO
	812	ENDDO
	813
[148]	814
[26]	815	END SUBROUTINE compute_gradiv
[15]	816
[151]	817	SUBROUTINE compute_divgrad(theta,divgrad_i,llb,lle)
[19]	818	USE icosa
[15]	819	IMPLICIT NONE
[151]	820	INTEGER,INTENT(IN) :: llb
	821	INTEGER,INTENT(IN) :: lle
	822	REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm)
	823	REAL(rstd),INTENT(OUT) :: divgrad_i(iim*jjm,llm)
	824	REAL(rstd) :: grad_e(3iimjjm,llb:lle)
[15]	825
[174]	826	INTEGER :: ij,l
[15]	827
[148]	828
[151]	829	DO l=llb,lle
[487]	830	!DIR$ SIMD
[174]	831	DO ij=ij_begin_ext,ij_end_ext
[15]	832
[174]	833	grad_e(ij+u_right,l)=-1/de(ij+u_right)(ne(ij,right)theta(ij,l)+ ne(ij+t_right,left)*theta(ij+t_right,l) )
[15]	834
[174]	835	grad_e(ij+u_lup,l)=-1/de(ij+u_lup)(ne(ij,lup)theta(ij,l)+ ne(ij+t_lup,rdown)*theta(ij+t_lup,l ))
[15]	836
[174]	837	grad_e(ij+u_ldown,l)=-1/de(ij+u_ldown)(ne(ij,ldown)theta(ij,l)+ne(ij+t_ldown,rup)*theta(ij+t_ldown,l) )
[15]	838
	839	ENDDO
	840	ENDDO
	841
	842
[151]	843	DO l=llb,lle
[487]	844	!DIR$ SIMD
[175]	845	DO ij=ij_begin,ij_end
[174]	846
	847	divgrad_i(ij,l)=1./Ai(ij)(ne(ij,right)grad_e(ij+u_right,l)*le(ij+u_right) + &
	848	ne(ij,rup)grad_e(ij+u_rup,l)le(ij+u_rup) + &
	849	ne(ij,lup)grad_e(ij+u_lup,l)le(ij+u_lup) + &
	850	ne(ij,left)grad_e(ij+u_left,l)le(ij+u_left) + &
	851	ne(ij,ldown)grad_e(ij+u_ldown,l)le(ij+u_ldown) + &
	852	ne(ij,rdown)grad_e(ij+u_rdown,l)le(ij+u_rdown))
[15]	853	ENDDO
	854	ENDDO
	855
[151]	856	DO l=llb,lle
[174]	857	DO ij=ij_begin,ij_end
	858	divgrad_i(ij,l)=-divgrad_i(ij,l)*cdivgrad
[15]	859	ENDDO
	860	ENDDO
	861
[26]	862	END SUBROUTINE compute_divgrad
[15]	863
	864
[151]	865	SUBROUTINE compute_gradrot(ue,gradrot_e,llb,lle)
[19]	866	USE icosa
[15]	867	IMPLICIT NONE
[151]	868	INTEGER,INTENT(IN) :: llb
	869	INTEGER,INTENT(IN) :: lle
	870	REAL(rstd),INTENT(IN) :: ue(iim3jjm,llm)
	871	REAL(rstd),INTENT(OUT) :: gradrot_e(iim3jjm,llm)
	872	REAL(rstd) :: rot_v(2iimjjm,llb:lle)
[15]	873
[174]	874	INTEGER :: ij,l
[15]	875
[151]	876	DO l=llb,lle
[487]	877	!DIR$ SIMD
[174]	878	DO ij=ij_begin_ext,ij_end_ext
[15]	879
[174]	880	rot_v(ij+z_up,l)= 1./Av(ij+z_up)( ne(ij,rup)ue(ij+u_rup,l)*de(ij+u_rup) &
	881	+ ne(ij+t_rup,left)ue(ij+t_rup+u_left,l)de(ij+t_rup+u_left) &
	882	- ne(ij,lup)ue(ij+u_lup,l)de(ij+u_lup) )
[15]	883
[174]	884	rot_v(ij+z_down,l) = 1./Av(ij+z_down)( ne(ij,ldown)ue(ij+u_ldown,l)*de(ij+u_ldown) &
	885	+ ne(ij+t_ldown,right)ue(ij+t_ldown+u_right,l)de(ij+t_ldown+u_right) &
	886	- ne(ij,rdown)ue(ij+u_rdown,l)de(ij+u_rdown) )
[15]	887
	888	ENDDO
	889	ENDDO
	890
[151]	891	DO l=llb,lle
[487]	892	!DIR$ SIMD
[174]	893	DO ij=ij_begin,ij_end
[15]	894
[174]	895	gradrot_e(ij+u_right,l)=1/le(ij+u_right)ne(ij,right)(rot_v(ij+z_rdown,l)-rot_v(ij+z_rup,l))
[15]	896
[174]	897	gradrot_e(ij+u_lup,l)=1/le(ij+u_lup)ne(ij,lup)(rot_v(ij+z_up,l)-rot_v(ij+z_lup,l))
[15]	898
[174]	899	gradrot_e(ij+u_ldown,l)=1/le(ij+u_ldown)ne(ij,ldown)(rot_v(ij+z_ldown,l)-rot_v(ij+z_down,l))
[15]	900
	901	ENDDO
	902	ENDDO
	903
[151]	904	DO l=llb,lle
[487]	905	!DIR$ SIMD
[174]	906	DO ij=ij_begin,ij_end
[15]	907
[174]	908	gradrot_e(ij+u_right,l)=-gradrot_e(ij+u_right,l)*cgradrot
	909	gradrot_e(ij+u_lup,l)=-gradrot_e(ij+u_lup,l)*cgradrot
	910	gradrot_e(ij+u_ldown,l)=-gradrot_e(ij+u_ldown,l)*cgradrot
[15]	911
	912	ENDDO
	913	ENDDO
	914
[26]	915	END SUBROUTINE compute_gradrot
[15]	916
	917
	918	END MODULE dissip_gcm_mod
	919

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

Download in other formats: