Context Navigation

source: codes/icosagcm/trunk/src/dissip/dissip_gcm.F90 @ 953

Last change on this file since 953 was 953, checked in by adurocher, 5 years ago
trunk : GPU implementation with OpenACC ( merge from glcp.idris.fr )
File size: 34.2 KB

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

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

Download in other formats: