Context Navigation

source: codes/icosagcm/trunk/src/dissip/dissip_gcm.f90 @ 648

Last change on this file since 648 was 648, checked in by aslmd, 7 years ago
added an option for simple Rayleigh friction at bottom (giant planets, but could be used for any case really)
File size: 27.2 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)
[524]	36	!$OMP THREADPRIVATE(rayleigh_shear)
[648]	37	REAL, SAVE :: rayleigh_tau, rayleigh_limlat
[524]	38	!$OMP THREADPRIVATE(rayleigh_tau)
[648]	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
[541]	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'
[648]	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
[648]	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
[541]	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
[541]	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)
[541]	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
[523]	451
	452	IF(rayleigh_friction_type>0) mintau=MIN(mintau,rayleigh_tau)
	453
[250]	454	IF(mintau>0) THEN
	455	itau_dissip=INT(mintau/dt)
	456	dtdissip=itau_dissip*dt
	457	ELSE
[295]	458	IF (is_master) PRINT *,"No dissipation time set, setting itau_dissip to 1000000000"
[250]	459	itau_dissip=100000000
	460	END IF
[148]	461	itau_dissip=MAX(1,itau_dissip)
[523]	462	IF (is_master) PRINT *,"rayleigh_tau",rayleigh_tau, "mintau ",mintau, &
	463	"itau_dissip",itau_dissip," dtdissip ",dtdissip
[15]	464
	465	END SUBROUTINE init_dissip
	466
	467
[523]	468	SUBROUTINE dissip(f_ps,f_mass,f_phis,f_geopot,f_theta_rhodz,f_ue, f_dtheta_rhodz,f_due)
[19]	469	USE icosa
[15]	470	USE theta2theta_rhodz_mod
[109]	471	USE pression_mod
	472	USE exner_mod
	473	USE geopotential_mod
[145]	474	USE trace
[148]	475	USE time_mod
[151]	476	USE omp_para
[15]	477	IMPLICIT NONE
[523]	478	TYPE(t_field),POINTER :: f_ps(:), f_mass(:), f_phis(:), f_geopot(:)
	479	TYPE(t_field),POINTER :: f_theta_rhodz(:), f_dtheta_rhodz(:)
	480	TYPE(t_field),POINTER :: f_ue(:), f_due(:)
[26]	481
[15]	482	REAL(rstd),POINTER :: due(:,:)
[109]	483	REAL(rstd),POINTER :: phi(:,:), ue(:,:)
[26]	484	REAL(rstd),POINTER :: due_diss1(:,:)
	485	REAL(rstd),POINTER :: due_diss2(:,:)
[387]	486	REAL(rstd),POINTER :: dtheta_rhodz(:,:,:)
[26]	487	REAL(rstd),POINTER :: dtheta_rhodz_diss(:,:)
[15]	488
[109]	489	INTEGER :: ind, shear
[648]	490	INTEGER :: l,ij,nn
[151]	491
	492	!$OMP BARRIER
[15]	493
[145]	494	CALL trace_start("dissip")
[26]	495	CALL gradiv(f_ue,f_due_diss1)
	496	CALL gradrot(f_ue,f_due_diss2)
[167]	497
	498	CALL divgrad_theta_rhodz(f_mass,f_theta_rhodz,f_dtheta_rhodz_diss)
[26]	499
[15]	500	DO ind=1,ndomain
[186]	501	IF (.NOT. assigned_domain(ind)) CYCLE
[15]	502	CALL swap_dimensions(ind)
	503	CALL swap_geometry(ind)
	504	due=f_due(ind)
[26]	505	due_diss1=f_due_diss1(ind)
	506	due_diss2=f_due_diss2(ind)
	507	dtheta_rhodz=f_dtheta_rhodz(ind)
	508	dtheta_rhodz_diss=f_dtheta_rhodz_diss(ind)
[151]	509
	510	DO l=ll_begin,ll_end
[487]	511	!DIR$ SIMD
[174]	512	DO ij=ij_begin,ij_end
[15]	513
[174]	514	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
	515	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
	516	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]	517
[387]	518	dtheta_rhodz(ij,l,1) = -0.5dtheta_rhodz_diss(ij,l)/tau_divgrad(l)itau_dissip
[26]	519	ENDDO
	520	ENDDO
[109]	521
[151]	522	! dtheta_rhodz=0
	523	! due=0
	524
[523]	525	IF(rayleigh_friction_type>0) THEN
[648]	526	IF(rayleigh_friction_type<99) THEN
[523]	527	phi=f_geopot(ind)
	528	ue=f_ue(ind)
[525]	529	DO l=ll_begin,ll_end
[523]	530	DO ij=ij_begin,ij_end
	531	CALL relax(t_right, u_right)
	532	CALL relax(t_lup, u_lup)
	533	CALL relax(t_ldown, u_ldown)
	534	ENDDO
	535	END DO
[648]	536	ELSE
	537	ue=f_ue(ind)
	538	DO ij=ij_begin,ij_end
	539	nn = ij+u_right
	540	IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN
	541	!print, "latitude", lat_e(nn)180./3.14159
	542	due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau)
	543	ENDIF
	544	nn = ij+u_lup
	545	IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN
	546	due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau)
	547	ENDIF
	548	nn = ij+u_ldown
	549	IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN
	550	due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau)
	551	ENDIF
	552	ENDDO
	553	ENDIF
[523]	554	END IF
[109]	555	END DO
	556
[145]	557	CALL trace_end("dissip")
	558
[151]	559	!$OMP BARRIER
	560
[109]	561	CONTAINS
[648]	562
[109]	563	SUBROUTINE relax(shift_t, shift_u)
	564	USE dcmip_initial_conditions_test_1_2_3
	565	REAL(rstd) :: z, ulon,ulat, lon,lat, & ! input to test2_schaer_mountain
	566	p,hyam,hybm,w,t,phis,ps,rho,q, & ! unused input/output to test2_schaer_mountain
	567	fz, u3d(3), uref
	568	REAL(rstd), PARAMETER :: zh=2e4,ztop=3e4 ! DCMIP values
	569	LOGICAL :: hybrid_eta
	570	INTEGER :: shift_u, shift_t, zcoords, nn
[174]	571	z = (phi(ij,l)+phi(ij+shift_t,l))/(2.*g)
[109]	572	IF(z>zh) THEN ! relax only in the sponge layer z>zh
[151]	573
[174]	574	nn = ij+shift_u
[109]	575	zcoords = 1 ; hybrid_eta = .FALSE. ! use z instead of p or hyam/hybm
[286]	576	CALL test2_schaer_mountain(lon_e(nn),lat_e(nn),p,z,zcoords,hybrid_eta, &
[109]	577	hyam,hybm,shear,ulon,ulat,w,t,phis,ps,rho,q)
	578	! u3d = ulonelon_e(nn,:) + ulatelat_e(nn,:)
	579	u3d = ulon*elon_e(nn,:) ! ulat=0
	580	uref = sum(u3d*ep_e(nn,:))
	581
	582	fz = sin((pi/2)*(z-zh)/(ztop-zh))
	583	fz = fz*fz/rayleigh_tau
[523]	584	due(nn,l) = due(nn,l) - itau_dissipfz(ue(nn,l)-uref)
	585	! fz = 1./rayleigh_tau
[109]	586	! due(nn,l) = due(nn,l) - fz*ue(nn,l)
	587	END IF
	588	END SUBROUTINE relax
[15]	589
	590	END SUBROUTINE dissip
	591
[26]	592	SUBROUTINE gradiv(f_ue,f_due)
[19]	593	USE icosa
[145]	594	USE trace
[151]	595	USE omp_para
[15]	596	IMPLICIT NONE
[26]	597	TYPE(t_field),POINTER :: f_ue(:)
	598	TYPE(t_field),POINTER :: f_due(:)
	599	REAL(rstd),POINTER :: ue(:,:)
	600	REAL(rstd),POINTER :: due(:,:)
	601	INTEGER :: ind
[174]	602	INTEGER :: it,l,ij
[26]	603
[145]	604	CALL trace_start("gradiv")
	605
[26]	606	DO ind=1,ndomain
[186]	607	IF (.NOT. assigned_domain(ind)) CYCLE
[26]	608	CALL swap_dimensions(ind)
	609	CALL swap_geometry(ind)
	610	ue=f_ue(ind)
	611	due=f_due(ind)
[151]	612	DO l = ll_begin, ll_end
[487]	613	!DIR$ SIMD
[174]	614	DO ij=ij_begin,ij_end
[151]	615	due(ij+u_right,l)=ue(ij+u_right,l)
	616	due(ij+u_lup,l)=ue(ij+u_lup,l)
	617	due(ij+u_ldown,l)=ue(ij+u_ldown,l)
	618	ENDDO
	619	ENDDO
[26]	620	ENDDO
[15]	621
[26]	622	DO it=1,nitergdiv
	623
[151]	624	CALL send_message(f_due,req_due)
	625	CALL wait_message(req_due)
[26]	626
	627	DO ind=1,ndomain
[186]	628	IF (.NOT. assigned_domain(ind)) CYCLE
[26]	629	CALL swap_dimensions(ind)
	630	CALL swap_geometry(ind)
	631	due=f_due(ind)
[151]	632	CALL compute_gradiv(due,due,ll_begin,ll_end)
[26]	633	ENDDO
	634	ENDDO
[15]	635
[145]	636	CALL trace_end("gradiv")
	637
[26]	638	END SUBROUTINE gradiv
	639
	640
	641	SUBROUTINE gradrot(f_ue,f_due)
	642	USE icosa
[145]	643	USE trace
[151]	644	USE omp_para
[26]	645	IMPLICIT NONE
	646	TYPE(t_field),POINTER :: f_ue(:)
	647	TYPE(t_field),POINTER :: f_due(:)
	648	REAL(rstd),POINTER :: ue(:,:)
	649	REAL(rstd),POINTER :: due(:,:)
[15]	650	INTEGER :: ind
[174]	651	INTEGER :: it,l,ij
[26]	652
[145]	653	CALL trace_start("gradrot")
	654
[26]	655	DO ind=1,ndomain
[186]	656	IF (.NOT. assigned_domain(ind)) CYCLE
[26]	657	CALL swap_dimensions(ind)
	658	CALL swap_geometry(ind)
	659	ue=f_ue(ind)
	660	due=f_due(ind)
[151]	661	DO l = ll_begin, ll_end
[487]	662	!DIR$ SIMD
[174]	663	DO ij=ij_begin,ij_end
[151]	664	due(ij+u_right,l)=ue(ij+u_right,l)
	665	due(ij+u_lup,l)=ue(ij+u_lup,l)
	666	due(ij+u_ldown,l)=ue(ij+u_ldown,l)
	667	ENDDO
	668	ENDDO
[26]	669	ENDDO
[15]	670
[26]	671	DO it=1,nitergrot
	672
[151]	673	CALL send_message(f_due,req_due)
	674	CALL wait_message(req_due)
[26]	675
	676	DO ind=1,ndomain
[186]	677	IF (.NOT. assigned_domain(ind)) CYCLE
[26]	678	CALL swap_dimensions(ind)
	679	CALL swap_geometry(ind)
	680	due=f_due(ind)
[151]	681	CALL compute_gradrot(due,due,ll_begin,ll_end)
[15]	682	ENDDO
	683
[26]	684	ENDDO
[15]	685
[145]	686	CALL trace_end("gradrot")
	687
[26]	688	END SUBROUTINE gradrot
	689
	690	SUBROUTINE divgrad(f_theta,f_dtheta)
	691	USE icosa
[145]	692	USE trace
[151]	693	USE omp_para
[26]	694	IMPLICIT NONE
	695	TYPE(t_field),POINTER :: f_theta(:)
	696	TYPE(t_field),POINTER :: f_dtheta(:)
	697	REAL(rstd),POINTER :: theta(:,:)
	698	REAL(rstd),POINTER :: dtheta(:,:)
	699	INTEGER :: ind
	700	INTEGER :: it
[145]	701
	702	CALL trace_start("divgrad")
[26]	703
	704	DO ind=1,ndomain
[186]	705	IF (.NOT. assigned_domain(ind)) CYCLE
[26]	706	CALL swap_dimensions(ind)
	707	CALL swap_geometry(ind)
	708	theta=f_theta(ind)
	709	dtheta=f_dtheta(ind)
	710	dtheta=theta
	711	ENDDO
[15]	712
[26]	713	DO it=1,niterdivgrad
	714
	715	CALL transfert_request(f_dtheta,req_i1)
	716
	717	DO ind=1,ndomain
[186]	718	IF (.NOT. assigned_domain(ind)) CYCLE
[26]	719	CALL swap_dimensions(ind)
	720	CALL swap_geometry(ind)
	721	dtheta=f_dtheta(ind)
[151]	722	CALL compute_divgrad(dtheta,dtheta,ll_begin,ll_end)
[26]	723	ENDDO
[15]	724
[26]	725	ENDDO
	726
[145]	727	CALL trace_end("divgrad")
	728
[26]	729	END SUBROUTINE divgrad
	730
[167]	731	SUBROUTINE divgrad_theta_rhodz(f_mass,f_theta_rhodz,f_dtheta_rhodz)
[151]	732	USE icosa
	733	USE trace
	734	USE omp_para
	735	IMPLICIT NONE
[167]	736	TYPE(t_field),POINTER :: f_mass(:)
[151]	737	TYPE(t_field),POINTER :: f_theta_rhodz(:)
	738	TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
	739
[167]	740	REAL(rstd),POINTER :: mass(:,:)
[387]	741	REAL(rstd),POINTER :: theta_rhodz(:,:,:)
[151]	742	REAL(rstd),POINTER :: dtheta_rhodz(:,:)
[26]	743
[151]	744	INTEGER :: ind
[174]	745	INTEGER :: it,l,ij
[151]	746
	747	CALL trace_start("divgrad")
[295]	748
[151]	749	DO ind=1,ndomain
[186]	750	IF (.NOT. assigned_domain(ind)) CYCLE
[151]	751	CALL swap_dimensions(ind)
	752	CALL swap_geometry(ind)
[167]	753	mass=f_mass(ind)
[151]	754	theta_rhodz=f_theta_rhodz(ind)
	755	dtheta_rhodz=f_dtheta_rhodz(ind)
	756	DO l = ll_begin, ll_end
[487]	757	!DIR$ SIMD
[174]	758	DO ij=ij_begin,ij_end
[387]	759	dtheta_rhodz(ij,l) = theta_rhodz(ij,l,1) / mass(ij,l)
[151]	760	ENDDO
	761	ENDDO
	762	ENDDO
	763
	764	DO it=1,niterdivgrad
	765
	766	CALL send_message(f_dtheta_rhodz,req_dtheta)
	767	CALL wait_message(req_dtheta)
	768	DO ind=1,ndomain
[186]	769	IF (.NOT. assigned_domain(ind)) CYCLE
[151]	770	CALL swap_dimensions(ind)
	771	CALL swap_geometry(ind)
	772	dtheta_rhodz=f_dtheta_rhodz(ind)
	773	CALL compute_divgrad(dtheta_rhodz,dtheta_rhodz,ll_begin,ll_end)
	774	ENDDO
	775
	776	ENDDO
	777
	778	DO ind=1,ndomain
[186]	779	IF (.NOT. assigned_domain(ind)) CYCLE
[151]	780	CALL swap_dimensions(ind)
	781	CALL swap_geometry(ind)
	782	dtheta_rhodz=f_dtheta_rhodz(ind)
[175]	783	mass=f_mass(ind)
[151]	784
	785	DO l = ll_begin, ll_end
[487]	786	!DIR$ SIMD
[174]	787	DO ij=ij_begin,ij_end
[167]	788	dtheta_rhodz(ij,l) = dtheta_rhodz(ij,l) * mass(ij,l)
[151]	789	ENDDO
	790	ENDDO
	791	ENDDO
	792
	793
	794	CALL trace_end("divgrad")
	795
	796	END SUBROUTINE divgrad_theta_rhodz
[15]	797
[151]	798	SUBROUTINE compute_gradiv(ue,gradivu_e,llb,lle)
[19]	799	USE icosa
[15]	800	IMPLICIT NONE
[151]	801	INTEGER,INTENT(IN) :: llb
	802	INTEGER,INTENT(IN) :: lle
	803	REAL(rstd),INTENT(IN) :: ue(iim3jjm,llm)
	804	REAL(rstd),INTENT(OUT) :: gradivu_e(iim3jjm,llm)
	805	REAL(rstd) :: divu_i(iim*jjm,llb:lle)
[15]	806
[174]	807	INTEGER :: ij,l
[15]	808
[151]	809	DO l=llb,lle
[487]	810	!DIR$ SIMD
[174]	811	DO ij=ij_begin,ij_end
	812	divu_i(ij,l)=1./Ai(ij)(ne(ij,right)ue(ij+u_right,l)*le(ij+u_right) + &
	813	ne(ij,rup)ue(ij+u_rup,l)le(ij+u_rup) + &
	814	ne(ij,lup)ue(ij+u_lup,l)le(ij+u_lup) + &
	815	ne(ij,left)ue(ij+u_left,l)le(ij+u_left) + &
	816	ne(ij,ldown)ue(ij+u_ldown,l)le(ij+u_ldown) + &
	817	ne(ij,rdown)ue(ij+u_rdown,l)le(ij+u_rdown))
[15]	818	ENDDO
	819	ENDDO
	820
[151]	821	DO l=llb,lle
[487]	822	!DIR$ SIMD
[174]	823	DO ij=ij_begin,ij_end
[15]	824
[174]	825	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]	826
[174]	827	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]	828
[174]	829	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]	830
	831	ENDDO
	832	ENDDO
	833
[151]	834	DO l=llb,lle
[487]	835	!DIR$ SIMD
[174]	836	DO ij=ij_begin,ij_end
	837	gradivu_e(ij+u_right,l)=-gradivu_e(ij+u_right,l)*cgraddiv
	838	gradivu_e(ij+u_lup,l)=-gradivu_e(ij+u_lup,l)*cgraddiv
	839	gradivu_e(ij+u_ldown,l)=-gradivu_e(ij+u_ldown,l)*cgraddiv
[15]	840	ENDDO
	841	ENDDO
	842
[148]	843
[26]	844	END SUBROUTINE compute_gradiv
[15]	845
[151]	846	SUBROUTINE compute_divgrad(theta,divgrad_i,llb,lle)
[19]	847	USE icosa
[15]	848	IMPLICIT NONE
[151]	849	INTEGER,INTENT(IN) :: llb
	850	INTEGER,INTENT(IN) :: lle
	851	REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm)
	852	REAL(rstd),INTENT(OUT) :: divgrad_i(iim*jjm,llm)
	853	REAL(rstd) :: grad_e(3iimjjm,llb:lle)
[15]	854
[174]	855	INTEGER :: ij,l
[15]	856
[148]	857
[151]	858	DO l=llb,lle
[487]	859	!DIR$ SIMD
[174]	860	DO ij=ij_begin_ext,ij_end_ext
[15]	861
[174]	862	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]	863
[174]	864	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]	865
[174]	866	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]	867
	868	ENDDO
	869	ENDDO
	870
	871
[151]	872	DO l=llb,lle
[487]	873	!DIR$ SIMD
[175]	874	DO ij=ij_begin,ij_end
[174]	875
	876	divgrad_i(ij,l)=1./Ai(ij)(ne(ij,right)grad_e(ij+u_right,l)*le(ij+u_right) + &
	877	ne(ij,rup)grad_e(ij+u_rup,l)le(ij+u_rup) + &
	878	ne(ij,lup)grad_e(ij+u_lup,l)le(ij+u_lup) + &
	879	ne(ij,left)grad_e(ij+u_left,l)le(ij+u_left) + &
	880	ne(ij,ldown)grad_e(ij+u_ldown,l)le(ij+u_ldown) + &
	881	ne(ij,rdown)grad_e(ij+u_rdown,l)le(ij+u_rdown))
[15]	882	ENDDO
	883	ENDDO
	884
[151]	885	DO l=llb,lle
[174]	886	DO ij=ij_begin,ij_end
	887	divgrad_i(ij,l)=-divgrad_i(ij,l)*cdivgrad
[15]	888	ENDDO
	889	ENDDO
	890
[26]	891	END SUBROUTINE compute_divgrad
[15]	892
	893
[151]	894	SUBROUTINE compute_gradrot(ue,gradrot_e,llb,lle)
[19]	895	USE icosa
[15]	896	IMPLICIT NONE
[151]	897	INTEGER,INTENT(IN) :: llb
	898	INTEGER,INTENT(IN) :: lle
	899	REAL(rstd),INTENT(IN) :: ue(iim3jjm,llm)
	900	REAL(rstd),INTENT(OUT) :: gradrot_e(iim3jjm,llm)
	901	REAL(rstd) :: rot_v(2iimjjm,llb:lle)
[15]	902
[174]	903	INTEGER :: ij,l
[15]	904
[151]	905	DO l=llb,lle
[487]	906	!DIR$ SIMD
[174]	907	DO ij=ij_begin_ext,ij_end_ext
[15]	908
[174]	909	rot_v(ij+z_up,l)= 1./Av(ij+z_up)( ne(ij,rup)ue(ij+u_rup,l)*de(ij+u_rup) &
	910	+ ne(ij+t_rup,left)ue(ij+t_rup+u_left,l)de(ij+t_rup+u_left) &
	911	- ne(ij,lup)ue(ij+u_lup,l)de(ij+u_lup) )
[15]	912
[174]	913	rot_v(ij+z_down,l) = 1./Av(ij+z_down)( ne(ij,ldown)ue(ij+u_ldown,l)*de(ij+u_ldown) &
	914	+ ne(ij+t_ldown,right)ue(ij+t_ldown+u_right,l)de(ij+t_ldown+u_right) &
	915	- ne(ij,rdown)ue(ij+u_rdown,l)de(ij+u_rdown) )
[15]	916
	917	ENDDO
	918	ENDDO
	919
[151]	920	DO l=llb,lle
[487]	921	!DIR$ SIMD
[174]	922	DO ij=ij_begin,ij_end
[15]	923
[174]	924	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]	925
[174]	926	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]	927
[174]	928	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]	929
	930	ENDDO
	931	ENDDO
	932
[151]	933	DO l=llb,lle
[487]	934	!DIR$ SIMD
[174]	935	DO ij=ij_begin,ij_end
[15]	936
[174]	937	gradrot_e(ij+u_right,l)=-gradrot_e(ij+u_right,l)*cgradrot
	938	gradrot_e(ij+u_lup,l)=-gradrot_e(ij+u_lup,l)*cgradrot
	939	gradrot_e(ij+u_ldown,l)=-gradrot_e(ij+u_ldown,l)*cgradrot
[15]	940
	941	ENDDO
	942	ENDDO
	943
[26]	944	END SUBROUTINE compute_gradrot
[15]	945
	946
	947	END MODULE dissip_gcm_mod
	948

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

Download in other formats: