Context Navigation

source: codes/icosagcm/trunk/src/dissip_gcm.f90 @ 59

Last change on this file since 59 was 59, checked in by ymipsl, 12 years ago

Correct non compliant operation in MPI_Allreduce

YM

File size: 21.5 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(:)
10	TYPE(t_field),POINTER,SAVE :: f_dtheta_diss(:)
11	TYPE(t_field),POINTER,SAVE :: f_dtheta_rhodz_diss(:)
12
13
14	INTEGER,SAVE :: nitergdiv=1
15	INTEGER,SAVE :: nitergrot=1
16	INTEGER,SAVE :: niterdivgrad=1
17
18	REAL,ALLOCATABLE,SAVE :: tau_graddiv(:)
19	REAL,ALLOCATABLE,SAVE :: tau_gradrot(:)
20	REAL,ALLOCATABLE,SAVE :: tau_divgrad(:)
21
22	REAL,SAVE :: cgraddiv
23	REAL,SAVE :: cgradrot
24	REAL,SAVE :: cdivgrad
25
26	INTEGER,SAVE :: idissip
27	REAL,SAVE :: dtdissip
28
29	PUBLIC init_dissip, dissip
30	CONTAINS
31
32	SUBROUTINE allocate_dissip
33	USE icosa
34	IMPLICIT NONE
35	CALL allocate_field(f_due_diss1,field_u,type_real,llm)
36	CALL allocate_field(f_due_diss2,field_u,type_real,llm)
37	CALL allocate_field(f_theta,field_t,type_real,llm)
38	CALL allocate_field(f_dtheta_diss,field_t,type_real,llm)
39	CALL allocate_field(f_dtheta_rhodz_diss,field_t,type_real,llm)
40
41	ALLOCATE(tau_graddiv(llm))
42	ALLOCATE(tau_gradrot(llm))
43	ALLOCATE(tau_divgrad(llm))
44	END SUBROUTINE allocate_dissip
45
46	SUBROUTINE init_dissip(dt)
47	USE icosa
48	USE disvert_mod
49	USE mpi_mod
50	USE mpipara
51
52	IMPLICIT NONE
53	REAL,INTENT(IN) :: dt
54
55	TYPE(t_field),POINTER :: f_u(:)
56	TYPE(t_field),POINTER :: f_du(:)
57	REAL(rstd),POINTER :: u(:)
58	REAL(rstd),POINTER :: du(:)
59	TYPE(t_field),POINTER :: f_theta(:)
60	TYPE(t_field),POINTER :: f_dtheta(:)
61	REAL(rstd),POINTER :: theta(:)
62	REAL(rstd),POINTER :: dtheta(:)
63	REAL(rstd) :: dumax,dumax1
64	REAL(rstd) :: dthetamax,dthetamax1
65	REAL(rstd) :: r
66	REAL(rstd) :: tau
67	REAL(rstd) :: zz, zvert, fact
68	INTEGER :: l
69
70
71	INTEGER :: i,j,n,ind,it,iter
72
73	CALL allocate_dissip
74	CALL allocate_field(f_u,field_u,type_real)
75	CALL allocate_field(f_du,field_u,type_real)
76	CALL allocate_field(f_theta,field_t,type_real)
77	CALL allocate_field(f_dtheta,field_t,type_real)
78
79	tau_graddiv(:)=5000
80	CALL getin("tau_graddiv",tau)
81	tau_graddiv(:)=tau/scale_factor
82
83	CALL getin("nitergdiv",nitergdiv)
84
85	tau_gradrot(:)=5000
86	CALL getin("tau_gradrot",tau_gradrot)
87	tau_gradrot(:)=tau/scale_factor
88
89	CALL getin("nitergrot",nitergrot)
90
91
92	tau_divgrad(:)=5000
93	CALL getin("tau_divgrad",tau)
94	tau_divgrad(:)=tau/scale_factor
95
96	CALL getin("niterdivgrad",niterdivgrad)
97
98	! CALL create_request(field_u,req_dissip)
99
100	! DO ind=1,ndomain
101	! DO i=ii_begin,ii_end
102	! CALL request_add_point(ind,i,jj_begin-1,req_dissip,rup)
103	! CALL request_add_point(ind,i+1,jj_begin-1,req_dissip,lup)
104	! ENDDO
105
106	! DO j=jj_begin,jj_end
107	! CALL request_add_point(ind,ii_end+1,j,req_dissip,left)
108	! CALL request_add_point(ind,ii_end+1,j-1,req_dissip,lup)
109	! ENDDO
110	!
111	! DO i=ii_begin,ii_end
112	! CALL request_add_point(ind,i,jj_end+1,req_dissip,ldown)
113	! CALL request_add_point(ind,i-1,jj_end+1,req_dissip,rdown)
114	! ENDDO
115
116	! DO j=jj_begin,jj_end
117	! CALL request_add_point(ind,ii_begin-1,j,req_dissip,right)
118	! CALL request_add_point(ind,ii_begin-1,j+1,req_dissip,rdown)
119	! ENDDO
120
121	! DO i=ii_begin+1,ii_end-1
122	! CALL request_add_point(ind,i,jj_begin,req_dissip,right)
123	! CALL request_add_point(ind,i,jj_end,req_dissip,right)
124	! ENDDO
125	!
126	! DO j=jj_begin+1,jj_end-1
127	! CALL request_add_point(ind,ii_begin,j,req_dissip,rup)
128	! CALL request_add_point(ind,ii_end,j,req_dissip,rup)
129	! ENDDO
130
131	! CALL request_add_point(ind,ii_begin+1,jj_begin,req_dissip,left)
132	! CALL request_add_point(ind,ii_begin,jj_begin+1,req_dissip,ldown)
133	! CALL request_add_point(ind,ii_begin+1,jj_end,req_dissip,left)
134	! CALL request_add_point(ind,ii_end,jj_begin+1,req_dissip,ldown)
135	!
136	! ENDDO
137
138
139	cgraddiv=-1
140	cdivgrad=-1
141	cgradrot=-1
142
143	CALL RANDOM_SEED()
144
145	DO ind=1,ndomain
146	CALL swap_dimensions(ind)
147	CALL swap_geometry(ind)
148	u=f_u(ind)
149
150	DO j=jj_begin,jj_end
151	DO i=ii_begin,ii_end
152	n=(j-1)*iim+i
153	CALL RANDOM_NUMBER(r)
154	u(n+u_right)=r-0.5
155	CALL RANDOM_NUMBER(r)
156	u(n+u_lup)=r-0.5
157	CALL RANDOM_NUMBER(r)
158	u(n+u_ldown)=r-0.5
159	ENDDO
160	ENDDO
161
162	ENDDO
163
164
165
166	DO it=1,20
167
168	dumax=0
169	DO iter=1,nitergdiv
170	CALL transfert_request(f_u,req_e1)
171	DO ind=1,ndomain
172	CALL swap_dimensions(ind)
173	CALL swap_geometry(ind)
174	u=f_u(ind)
175	du=f_du(ind)
176	CALL compute_gradiv(u,du,1)
177	u=du
178	ENDDO
179	ENDDO
180
181	CALL transfert_request(f_du,req_e1)
182
183	DO ind=1,ndomain
184	CALL swap_dimensions(ind)
185	CALL swap_geometry(ind)
186	u=f_u(ind)
187	du=f_du(ind)
188
189	DO j=jj_begin,jj_end
190	DO i=ii_begin,ii_end
191	n=(j-1)*iim+i
192	if (le(n+u_right)>1e-100) dumax=MAX(dumax,ABS(du(n+u_right)))
193	if (le(n+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(n+u_lup)))
194	if (le(n+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(n+u_ldown)))
195	ENDDO
196	ENDDO
197	ENDDO
198
199	IF (using_mpi) THEN
200	CALL MPI_ALLREDUCE(dumax,dumax1,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
201	dumax=dumax1
202	ENDIF
203
204	DO ind=1,ndomain
205	CALL swap_dimensions(ind)
206	CALL swap_geometry(ind)
207	u=f_u(ind)
208	du=f_du(ind)
209	u=du/dumax
210	ENDDO
211	PRINT *,"gradiv : it :",it ,": dumax",dumax
212
213	ENDDO
214	PRINT *,"gradiv : dumax",dumax
215
216	cgraddiv=dumax**(-1./nitergdiv)
217	PRINT *,"cgraddiv : ",cgraddiv
218
219	DO ind=1,ndomain
220	CALL swap_dimensions(ind)
221	CALL swap_geometry(ind)
222	u=f_u(ind)
223
224	DO j=jj_begin,jj_end
225	DO i=ii_begin,ii_end
226	n=(j-1)*iim+i
227	CALL RANDOM_NUMBER(r)
228	u(n+u_right)=r-0.5
229	CALL RANDOM_NUMBER(r)
230	u(n+u_lup)=r-0.5
231	CALL RANDOM_NUMBER(r)
232	u(n+u_ldown)=r-0.5
233	ENDDO
234	ENDDO
235
236	ENDDO
237
238
239	DO it=1,20
240
241	dumax=0
242	DO iter=1,nitergrot
243	CALL transfert_request(f_u,req_e1)
244	DO ind=1,ndomain
245	CALL swap_dimensions(ind)
246	CALL swap_geometry(ind)
247	u=f_u(ind)
248	du=f_du(ind)
249	CALL compute_gradrot(u,du,1)
250	u=du
251	ENDDO
252	ENDDO
253
254	CALL transfert_request(f_du,req_e1)
255
256	DO ind=1,ndomain
257	CALL swap_dimensions(ind)
258	CALL swap_geometry(ind)
259	u=f_u(ind)
260	du=f_du(ind)
261
262	DO j=jj_begin,jj_end
263	DO i=ii_begin,ii_end
264	n=(j-1)*iim+i
265	if (le(n+u_right)>1e-100) dumax=MAX(dumax,ABS(du(n+u_right)))
266	if (le(n+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(n+u_lup)))
267	if (le(n+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(n+u_ldown)))
268	ENDDO
269	ENDDO
270	ENDDO
271
272	IF (using_mpi) THEN
273	CALL MPI_ALLREDUCE(dumax,dumax1,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
274	dumax=dumax1
275	ENDIF
276
277	DO ind=1,ndomain
278	CALL swap_dimensions(ind)
279	CALL swap_geometry(ind)
280	u=f_u(ind)
281	du=f_du(ind)
282	u=du/dumax
283	ENDDO
284
285	PRINT *,"gradrot : it :",it ,": dumax",dumax
286
287	ENDDO
288	PRINT *,"gradrot : dumax",dumax
289
290	cgradrot=dumax**(-1./nitergrot)
291	PRINT *,"cgradrot : ",cgradrot
292
293
294
295	DO ind=1,ndomain
296	CALL swap_dimensions(ind)
297	CALL swap_geometry(ind)
298	theta=f_theta(ind)
299
300	DO j=jj_begin,jj_end
301	DO i=ii_begin,ii_end
302	n=(j-1)*iim+i
303	CALL RANDOM_NUMBER(r)
304	theta(n)=r-0.5
305	ENDDO
306	ENDDO
307
308	ENDDO
309
310	DO it=1,20
311
312	dthetamax=0
313	DO iter=1,niterdivgrad
314	CALL transfert_request(f_theta,req_i1)
315	DO ind=1,ndomain
316	CALL swap_dimensions(ind)
317	CALL swap_geometry(ind)
318	theta=f_theta(ind)
319	dtheta=f_dtheta(ind)
320	CALL compute_divgrad(theta,dtheta,1)
321	theta=dtheta
322	ENDDO
323	ENDDO
324	! CALL writefield("divgrad",f_dtheta)
325
326	CALL transfert_request(f_dtheta,req_i1)
327
328	DO ind=1,ndomain
329	CALL swap_dimensions(ind)
330	CALL swap_geometry(ind)
331	theta=f_theta(ind)
332	dtheta=f_dtheta(ind)
333
334	DO j=jj_begin,jj_end
335	DO i=ii_begin,ii_end
336	n=(j-1)*iim+i
337	dthetamax=MAX(dthetamax,ABS(dtheta(n)))
338	ENDDO
339	ENDDO
340	ENDDO
341	IF (using_mpi) THEN
342	CALL MPI_ALLREDUCE(dthetamax,dthetamax1,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
343	dthetamax=dthetamax1
344	ENDIF
345	PRINT *,"divgrad : it :",it ,": dthetamax",dthetamax
346
347	DO ind=1,ndomain
348	CALL swap_dimensions(ind)
349	CALL swap_geometry(ind)
350	theta=f_theta(ind)
351	dtheta=f_dtheta(ind)
352	theta=dtheta/dthetamax
353	ENDDO
354	ENDDO
355
356	! CALL writefield("divgrad",f_dtheta)
357	PRINT *,"divgrad : divgrad",dthetamax
358
359	cdivgrad=dthetamax**(-1./niterdivgrad)
360	PRINT *,"cdivgrad : ",cdivgrad
361
362
363
364
365
366
367
368	fact=2
369	DO l=1,llm
370	zz= 1. - preff/presnivs(l)
371	zvert=fact-(fact-1)/(1+zz*zz)
372	tau_graddiv(l) = tau_graddiv(l)/zvert
373	tau_gradrot(l) = tau_gradrot(l)/zvert
374	tau_divgrad(l) = tau_divgrad(l)/zvert
375	ENDDO
376
377	! idissip=INT(MIN(tetagdiv,tetagrot)/(2.*dt))
378	! idissip=MAX(1,idissip)
379	! dtdissip=idissip*dt
380	! PRINT *,"idissip",idissip," dtdissip ",dtdissip
381
382	END SUBROUTINE init_dissip
383
384
385	SUBROUTINE dissip(f_ue,f_due,f_ps,f_theta_rhodz,f_dtheta_rhodz)
386	USE icosa
387	USE theta2theta_rhodz_mod
388	IMPLICIT NONE
389	TYPE(t_field),POINTER :: f_ue(:)
390	TYPE(t_field),POINTER :: f_due(:)
391	TYPE(t_field),POINTER :: f_ps(:)
392	TYPE(t_field),POINTER :: f_theta_rhodz(:)
393	TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
394
395	REAL(rstd),POINTER :: due(:,:)
396	REAL(rstd),POINTER :: due_diss1(:,:)
397	REAL(rstd),POINTER :: due_diss2(:,:)
398	REAL(rstd),POINTER :: dtheta_rhodz(:,:)
399	REAL(rstd),POINTER :: dtheta_rhodz_diss(:,:)
400
401	INTEGER :: ind
402	INTEGER :: l,i,j,n
403
404	CALL gradiv(f_ue,f_due_diss1)
405	CALL gradrot(f_ue,f_due_diss2)
406	CALL theta_rhodz2theta(f_ps,f_theta_rhodz,f_theta)
407	CALL divgrad(f_theta,f_dtheta_diss)
408	CALL theta2theta_rhodz(f_ps,f_dtheta_diss,f_dtheta_rhodz_diss)
409
410	DO ind=1,ndomain
411	CALL swap_dimensions(ind)
412	CALL swap_geometry(ind)
413	due=f_due(ind)
414	due_diss1=f_due_diss1(ind)
415	due_diss2=f_due_diss2(ind)
416	dtheta_rhodz=f_dtheta_rhodz(ind)
417	dtheta_rhodz_diss=f_dtheta_rhodz_diss(ind)
418
419	DO l=1,llm
420	DO j=jj_begin,jj_end
421	DO i=ii_begin,ii_end
422	n=(j-1)*iim+i
423
424	due(n+u_right,l)=due(n+u_right,l)-0.5*(due_diss1(n+u_right,l)/tau_graddiv(l)+due_diss2(n+u_right,l)/tau_gradrot(l))
425	due(n+u_lup,l)=due(n+u_lup,l) -0.5*(due_diss1(n+u_lup,l) /tau_graddiv(l)+due_diss2(n+u_lup,l) /tau_gradrot(l))
426	due(n+u_ldown,l)=due(n+u_ldown,l)-0.5*(due_diss1(n+u_ldown,l)/tau_graddiv(l)+due_diss2(n+u_ldown,l)/tau_gradrot(l))
427
428	dtheta_rhodz(n,l)=dtheta_rhodz(n,l)-0.5*dtheta_rhodz_diss(n,l)/tau_divgrad(l)
429
430	ENDDO
431	ENDDO
432	ENDDO
433	ENDDO
434
435	END SUBROUTINE dissip
436
437	SUBROUTINE gradiv(f_ue,f_due)
438	USE icosa
439	IMPLICIT NONE
440	TYPE(t_field),POINTER :: f_ue(:)
441	TYPE(t_field),POINTER :: f_due(:)
442	REAL(rstd),POINTER :: ue(:,:)
443	REAL(rstd),POINTER :: due(:,:)
444	INTEGER :: ind
445	INTEGER :: it
446
447	DO ind=1,ndomain
448	CALL swap_dimensions(ind)
449	CALL swap_geometry(ind)
450	ue=f_ue(ind)
451	due=f_due(ind)
452	due=ue
453	ENDDO
454
455	DO it=1,nitergdiv
456
457	CALL transfert_request(f_due,req_e1)
458
459	DO ind=1,ndomain
460	CALL swap_dimensions(ind)
461	CALL swap_geometry(ind)
462	due=f_due(ind)
463	CALL compute_gradiv(due,due,llm)
464	ENDDO
465	ENDDO
466
467	END SUBROUTINE gradiv
468
469
470	SUBROUTINE gradrot(f_ue,f_due)
471	USE icosa
472	IMPLICIT NONE
473	TYPE(t_field),POINTER :: f_ue(:)
474	TYPE(t_field),POINTER :: f_due(:)
475	REAL(rstd),POINTER :: ue(:,:)
476	REAL(rstd),POINTER :: due(:,:)
477	INTEGER :: ind
478	INTEGER :: it
479
480	DO ind=1,ndomain
481	CALL swap_dimensions(ind)
482	CALL swap_geometry(ind)
483	ue=f_ue(ind)
484	due=f_due(ind)
485	due=ue
486	ENDDO
487
488	DO it=1,nitergrot
489
490	CALL transfert_request(f_due,req_e1)
491
492	DO ind=1,ndomain
493	CALL swap_dimensions(ind)
494	CALL swap_geometry(ind)
495	due=f_due(ind)
496	CALL compute_gradrot(due,due,llm)
497	ENDDO
498
499	ENDDO
500
501	END SUBROUTINE gradrot
502
503	SUBROUTINE divgrad(f_theta,f_dtheta)
504	USE icosa
505	IMPLICIT NONE
506	TYPE(t_field),POINTER :: f_theta(:)
507	TYPE(t_field),POINTER :: f_dtheta(:)
508	REAL(rstd),POINTER :: theta(:,:)
509	REAL(rstd),POINTER :: dtheta(:,:)
510	INTEGER :: ind
511	INTEGER :: it
512
513	DO ind=1,ndomain
514	CALL swap_dimensions(ind)
515	CALL swap_geometry(ind)
516	theta=f_theta(ind)
517	dtheta=f_dtheta(ind)
518	dtheta=theta
519	ENDDO
520
521	DO it=1,niterdivgrad
522
523	CALL transfert_request(f_dtheta,req_i1)
524
525	DO ind=1,ndomain
526	CALL swap_dimensions(ind)
527	CALL swap_geometry(ind)
528	dtheta=f_dtheta(ind)
529	CALL compute_divgrad(dtheta,dtheta,llm)
530	ENDDO
531
532	ENDDO
533
534	END SUBROUTINE divgrad
535
536
537
538	! SUBROUTINE compute_dissip(ue,due,ps,theta_rhodz,dtheta_rhodz)
539	! USE icosa
540	! USE theta2theta_rhodz_mod
541	! IMPLICIT NONE
542	! REAL(rstd) :: ue(3iimjjm,llm)
543	! REAL(rstd) :: due(3iimjjm,llm)
544	! REAL(rstd) :: ps(iim*jjm)
545	! REAL(rstd) :: theta_rhodz(iim*jjm,llm)
546	! REAL(rstd) :: dtheta_rhodz(iim*jjm,llm)
547	!
548	! REAL(rstd),SAVE,ALLOCATABLE :: theta(:,:)
549	! REAL(rstd),SAVE,ALLOCATABLE :: du_dissip(:,:)
550	! REAL(rstd),SAVE,ALLOCATABLE :: dtheta_dissip(:,:)
551	! REAL(rstd),SAVE,ALLOCATABLE :: dtheta_rhodz_dissip(:,:)
552	!
553	! INTEGER :: ind
554	! INTEGER :: l,i,j,n
555	!
556	!!$OMP BARRIER
557	!!$OMP MASTER
558	! ALLOCATE(theta(iim*jjm,llm))
559	! ALLOCATE(du_dissip(3iimjjm,llm))
560	! ALLOCATE(dtheta_dissip(iim*jjm,llm))
561	! ALLOCATE(dtheta_rhodz_dissip(iim*jjm,llm))
562	!!$OMP END MASTER
563	!!$OMP BARRIER
564	!
565	! CALL gradiv(ue,du_dissip,llm)
566	! DO l=1,llm
567	!!$OMP DO
568	! DO j=jj_begin,jj_end
569	! DO i=ii_begin,ii_end
570	! n=(j-1)*iim+i
571	! due(n+u_right,l)=due(n+u_right,l)+du_dissip(n+u_right,l)/tau_graddiv(l)*0.5
572	! due(n+u_lup,l)=due(n+u_lup,l)+du_dissip(n+u_lup,l)/tau_graddiv(l)*0.5
573	! due(n+u_ldown,l)=due(n+u_ldown,l)+du_dissip(n+u_ldown,l)/tau_graddiv(l)*0.5
574	! ENDDO
575	! ENDDO
576	! ENDDO
577	!
578	! CALL gradrot(ue,du_dissip,llm)
579	!
580	! DO l=1,llm
581	!!$OMP DO
582	! DO j=jj_begin,jj_end
583	! DO i=ii_begin,ii_end
584	! n=(j-1)*iim+i
585	! due(n+u_right,l)=due(n+u_right,l)+du_dissip(n+u_right,l)/tau_gradrot(l)*0.5
586	! due(n+u_lup,l)=due(n+u_lup,l)+du_dissip(n+u_lup,l)/tau_gradrot(l)*0.5
587	! due(n+u_ldown,l)=due(n+u_ldown,l)+du_dissip(n+u_ldown,l)/tau_gradrot(l)*0.5
588	! ENDDO
589	! ENDDO
590	! ENDDO
591	!
592	! CALL compute_theta_rhodz2theta(ps,theta_rhodz,theta,1)
593	! CALL divgrad(theta,dtheta_dissip,llm)
594	! CALL compute_theta2theta_rhodz(ps,dtheta_dissip,dtheta_rhodz_dissip,0)
595	!
596	! DO l=1,llm
597	!!$OMP DO
598	! DO j=jj_begin,jj_end
599	! DO i=ii_begin,ii_end
600	! n=(j-1)*iim+i
601	! dtheta_rhodz(n,l)=dtheta_rhodz(n,l)+dtheta_rhodz_dissip(n,l)/tau_divgrad(l)*0.5
602	! ENDDO
603	! ENDDO
604	! ENDDO
605	!
606	!!$OMP BARRIER
607	!!$OMP MASTER
608	! DEALLOCATE(theta)
609	! DEALLOCATE(du_dissip)
610	! DEALLOCATE(dtheta_dissip)
611	! DEALLOCATE(dtheta_rhodz_dissip)
612	!!$OMP END MASTER
613	!!$OMP BARRIER
614	!
615	! END SUBROUTINE compute_dissip
616
617
618	SUBROUTINE compute_gradiv(ue,gradivu_e,ll)
619	USE icosa
620	IMPLICIT NONE
621	INTEGER,INTENT(IN) :: ll
622	REAL(rstd),INTENT(IN) :: ue(iim3jjm,ll)
623	REAL(rstd),INTENT(OUT) :: gradivu_e(iim3jjm,ll)
624	REAL(rstd),SAVE,ALLOCATABLE :: divu_i(:,:)
625
626	INTEGER :: i,j,n,l
627
628	!$OMP BARRIER
629	!$OMP MASTER
630	ALLOCATE(divu_i(iim*jjm,ll))
631	!$OMP END MASTER
632	!$OMP BARRIER
633
634	DO l=1,ll
635	!$OMP DO
636	DO j=jj_begin,jj_end
637	DO i=ii_begin,ii_end
638	n=(j-1)*iim+i
639	divu_i(n,l)=1./Ai(n)(ne(n,right)ue(n+u_right,l)*le(n+u_right) + &
640	ne(n,rup)ue(n+u_rup,l)le(n+u_rup) + &
641	ne(n,lup)ue(n+u_lup,l)le(n+u_lup) + &
642	ne(n,left)ue(n+u_left,l)le(n+u_left) + &
643	ne(n,ldown)ue(n+u_ldown,l)le(n+u_ldown) + &
644	ne(n,rdown)ue(n+u_rdown,l)le(n+u_rdown))
645	ENDDO
646	ENDDO
647	ENDDO
648
649	DO l=1,ll
650	!$OMP DO
651	DO j=jj_begin,jj_end
652	DO i=ii_begin,ii_end
653
654	n=(j-1)*iim+i
655
656	gradivu_e(n+u_right,l)=-1/de(n+u_right)(ne(n,right)divu_i(n,l)+ ne(n+t_right,left)*divu_i(n+t_right,l) )
657
658	gradivu_e(n+u_lup,l)=-1/de(n+u_lup)(ne(n,lup)divu_i(n,l)+ ne(n+t_lup,rdown)*divu_i(n+t_lup,l))
659
660	gradivu_e(n+u_ldown,l)=-1/de(n+u_ldown)(ne(n,ldown)divu_i(n,l)+ne(n+t_ldown,rup)*divu_i(n+t_ldown,l) )
661
662	ENDDO
663	ENDDO
664	ENDDO
665
666	DO l=1,ll
667	!$OMP DO
668	DO j=jj_begin,jj_end
669	DO i=ii_begin,ii_end
670	n=(j-1)*iim+i
671	gradivu_e(n+u_right,l)=-gradivu_e(n+u_right,l)*cgraddiv
672	gradivu_e(n+u_lup,l)=-gradivu_e(n+u_lup,l)*cgraddiv
673	gradivu_e(n+u_ldown,l)=-gradivu_e(n+u_ldown,l)*cgraddiv
674	ENDDO
675	ENDDO
676	ENDDO
677
678	!$OMP BARRIER
679	!$OMP MASTER
680	DEALLOCATE(divu_i)
681	!$OMP END MASTER
682	!$OMP BARRIER
683
684
685	END SUBROUTINE compute_gradiv
686
687	SUBROUTINE compute_divgrad(theta,divgrad_i,ll)
688	USE icosa
689	IMPLICIT NONE
690	INTEGER,INTENT(IN) :: ll
691	REAL(rstd),INTENT(IN) :: theta(iim*jjm,ll)
692	REAL(rstd),INTENT(OUT) :: divgrad_i(iim*jjm,ll)
693	REAL(rstd),SAVE,ALLOCATABLE :: grad_e(:,:)
694
695	INTEGER :: i,j,n,l
696
697	!$OMP BARRIER
698	!$OMP MASTER
699	ALLOCATE(grad_e(3iimjjm,ll))
700	!$OMP END MASTER
701	!$OMP BARRIER
702
703	DO l=1,ll
704	!$OMP DO
705	DO j=jj_begin-1,jj_end+1
706	DO i=ii_begin-1,ii_end+1
707
708	n=(j-1)*iim+i
709
710	grad_e(n+u_right,l)=-1/de(n+u_right)(ne(n,right)theta(n,l)+ ne(n+t_right,left)*theta(n+t_right,l) )
711
712	grad_e(n+u_lup,l)=-1/de(n+u_lup)(ne(n,lup)theta(n,l)+ ne(n+t_lup,rdown)*theta(n+t_lup,l ))
713
714	grad_e(n+u_ldown,l)=-1/de(n+u_ldown)(ne(n,ldown)theta(n,l)+ne(n+t_ldown,rup)*theta(n+t_ldown,l) )
715
716	ENDDO
717	ENDDO
718	ENDDO
719
720
721	DO l=1,ll
722	!$OMP DO
723	DO j=jj_begin,jj_end
724	DO i=ii_begin,ii_end
725	n=(j-1)*iim+i
726	divgrad_i(n,l)=1./Ai(n)(ne(n,right)grad_e(n+u_right,l)*le(n+u_right) + &
727	ne(n,rup)grad_e(n+u_rup,l)le(n+u_rup) + &
728	ne(n,lup)grad_e(n+u_lup,l)le(n+u_lup) + &
729	ne(n,left)grad_e(n+u_left,l)le(n+u_left) + &
730	ne(n,ldown)grad_e(n+u_ldown,l)le(n+u_ldown) + &
731	ne(n,rdown)grad_e(n+u_rdown,l)le(n+u_rdown))
732	ENDDO
733	ENDDO
734	ENDDO
735
736	DO l=1,ll
737	!$OMP DO
738	DO j=jj_begin,jj_end
739	DO i=ii_begin,ii_end
740	n=(j-1)*iim+i
741	divgrad_i(n,l)=-divgrad_i(n,l)*cdivgrad
742	ENDDO
743	ENDDO
744	ENDDO
745
746	!$OMP BARRIER
747	!$OMP MASTER
748	DEALLOCATE(grad_e)
749	!$OMP END MASTER
750	!$OMP BARRIER
751
752	END SUBROUTINE compute_divgrad
753
754
755	SUBROUTINE compute_gradrot(ue,gradrot_e,ll)
756	USE icosa
757	IMPLICIT NONE
758	INTEGER,INTENT(IN) :: ll
759	REAL(rstd),INTENT(IN) :: ue(iim3jjm,ll)
760	REAL(rstd),INTENT(OUT) :: gradrot_e(iim3jjm,ll)
761	REAL(rstd),SAVE,ALLOCATABLE :: rot_v(:,:)
762
763	INTEGER :: i,j,n,l
764
765	!$OMP BARRIER
766	!$OMP MASTER
767	ALLOCATE(rot_v(2iimjjm,ll))
768	!$OMP END MASTER
769	!$OMP BARRIER
770
771	DO l=1,ll
772	!$OMP DO
773	DO j=jj_begin-1,jj_end+1
774	DO i=ii_begin-1,ii_end+1
775	n=(j-1)*iim+i
776
777	rot_v(n+z_up,l)= 1./Av(n+z_up)( ne(n,rup)ue(n+u_rup,l)*de(n+u_rup) &
778	+ ne(n+t_rup,left)ue(n+t_rup+u_left,l)de(n+t_rup+u_left) &
779	- ne(n,lup)ue(n+u_lup,l)de(n+u_lup) )
780
781	rot_v(n+z_down,l) = 1./Av(n+z_down)( ne(n,ldown)ue(n+u_ldown,l)*de(n+u_ldown) &
782	+ ne(n+t_ldown,right)ue(n+t_ldown+u_right,l)de(n+t_ldown+u_right) &
783	- ne(n,rdown)ue(n+u_rdown,l)de(n+u_rdown) )
784
785	ENDDO
786	ENDDO
787	ENDDO
788
789	DO l=1,ll
790	!$OMP DO
791	DO j=jj_begin,jj_end
792	DO i=ii_begin,ii_end
793	n=(j-1)*iim+i
794
795	gradrot_e(n+u_right,l)=1/le(n+u_right)ne(n,right)(rot_v(n+z_rdown,l)-rot_v(n+z_rup,l))
796
797	gradrot_e(n+u_lup,l)=1/le(n+u_lup)ne(n,lup)(rot_v(n+z_up,l)-rot_v(n+z_lup,l))
798
799	gradrot_e(n+u_ldown,l)=1/le(n+u_ldown)ne(n,ldown)(rot_v(n+z_ldown,l)-rot_v(n+z_down,l))
800
801	ENDDO
802	ENDDO
803	ENDDO
804
805	DO l=1,ll
806	!$OMP DO
807	DO j=jj_begin,jj_end
808	DO i=ii_begin,ii_end
809	n=(j-1)*iim+i
810
811	gradrot_e(n+u_right,l)=-gradrot_e(n+u_right,l)*cgradrot
812	gradrot_e(n+u_lup,l)=-gradrot_e(n+u_lup,l)*cgradrot
813	gradrot_e(n+u_ldown,l)=-gradrot_e(n+u_ldown,l)*cgradrot
814
815	ENDDO
816	ENDDO
817	ENDDO
818
819	!$OMP BARRIER
820	!$OMP MASTER
821	DEALLOCATE(rot_v)
822	!$OMP END MASTER
823	!$OMP BARRIER
824
825	END SUBROUTINE compute_gradrot
826
827
828	END MODULE dissip_gcm_mod
829

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

Download in other formats: