Context Navigation

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

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

Put time variable : dt, itaumax, write_period, itau_out in the time module

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

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

Download in other formats: