Context Navigation

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

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

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

Download in other formats: