Context Navigation

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

Last change on this file since 915 was 856, checked in by dubos, 5 years ago
devel : towards Fortran driver for unstructured/LAM meshes
File size: 26.7 KB

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

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

Download in other formats: