Context Navigation

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

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

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

Download in other formats: