Context Navigation

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

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

Bug fixe for output tracer

YM

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

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

Download in other formats: