Context Navigation

dissip_gcm.f90 @ 3843

Last change on this file since 3843 was 3810, checked in by ymipsl, 11 years ago
Add DYNAMICO in aquaplanet configuration YM
File size: 26.0 KB

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

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

Download in other formats:

Original Format