Context Navigation

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

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

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

Download in other formats: