Context Navigation

source: codes/icosagcm/trunk/src/caldyn_gcm.f90 @ 12

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

dynamico tree creation

YM

File size: 22.8 KB

Line
1	MODULE caldyn_gcm_mod
2	USE genmod
3	USE field_mod
4	USE transfert_mod
5
6	TYPE(t_field),POINTER :: f_out(:)
7	REAL(rstd),POINTER :: out(:,:)
8	TYPE(t_field),POINTER :: f_out_u(:)
9	REAL(rstd),POINTER :: out_u(:,:)
10	TYPE(t_field),POINTER :: f_out_z(:)
11	REAL(rstd),POINTER :: out_z(:,:)
12
13	CONTAINS
14
15	SUBROUTINE allocate_caldyn
16	USE domain_mod
17	USE dimensions
18	USE geometry
19	USE metric
20	IMPLICIT NONE
21	INTEGER :: ind,i,j
22
23	CALL allocate_field(f_out,field_t,type_real,llm)
24	CALL allocate_field(f_out_u,field_u,type_real,llm)
25	CALL allocate_field(f_out_z,field_z,type_real,llm)
26
27	END SUBROUTINE allocate_caldyn
28
29	SUBROUTINE swap_caldyn(ind)
30	IMPLICIT NONE
31	INTEGER,INTENT(IN) :: ind
32	out=f_out(ind)
33	out_u=f_out_u(ind)
34	out_z=f_out_z(ind)
35
36	END SUBROUTINE swap_caldyn
37
38	SUBROUTINE check_mass_conservation(f_ps,f_dps)
39	USE domain_mod
40	USE dimensions
41	USE geometry
42	IMPLICIT NONE
43	TYPE(t_field),POINTER :: f_ps(:)
44	TYPE(t_field),POINTER :: f_dps(:)
45	REAL(rstd),POINTER :: ps(:)
46	REAL(rstd),POINTER :: dps(:)
47	REAL(rstd) :: mass_tot,dmass_tot
48	INTEGER :: ind,i,j,ij
49
50	mass_tot=0
51	dmass_tot=0
52
53	CALL transfert_request(f_dps,req_i1)
54	CALL transfert_request(f_ps,req_i1)
55
56	DO ind=1,ndomain
57	CALL swap_dimensions(ind)
58	CALL swap_geometry(ind)
59
60	ps=f_ps(ind)
61	dps=f_dps(ind)
62
63	DO j=jj_begin,jj_end
64	DO i=ii_begin,ii_end
65	ij=(j-1)*iim+i
66	IF (domain(ind)%own(i,j)) THEN
67	mass_tot=mass_tot+ps(ij)*Ai(ij)/g
68	dmass_tot=dmass_tot+dps(ij)*Ai(ij)/g
69	ENDIF
70	ENDDO
71	ENDDO
72
73	ENDDO
74	PRINT*, "mass_tot ", mass_tot," dmass_tot ",dmass_tot
75
76	END SUBROUTINE check_mass_conservation
77
78
79
80	SUBROUTINE caldyn(f_phis, f_ps, f_theta_rhodz, f_u, f_dps, f_dtheta_rhodz, f_du)
81	USE domain_mod
82	USE dimensions
83	USE grid_param
84	USE geometry
85	USE metric
86	USE write_field
87	USE vorticity_mod
88	USE kinetic_mod
89	IMPLICIT NONE
90	TYPE(t_field),POINTER :: f_phis(:)
91	TYPE(t_field),POINTER :: f_ps(:)
92	TYPE(t_field),POINTER :: f_theta_rhodz(:)
93	TYPE(t_field),POINTER :: f_u(:)
94	TYPE(t_field),POINTER :: f_dps(:)
95	TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
96	TYPE(t_field),POINTER :: f_du(:)
97
98	REAL(rstd),POINTER :: phis(:)
99	REAL(rstd),POINTER :: ps(:)
100	REAL(rstd),POINTER :: theta_rhodz(:,:)
101	REAL(rstd),POINTER :: u(:,:)
102	REAL(rstd),POINTER :: dps(:)
103	REAL(rstd),POINTER :: dtheta_rhodz(:,:)
104	REAL(rstd),POINTER :: du(:,:)
105	INTEGER :: ind
106	INTEGER,SAVE :: it=0
107
108	CALL transfert_request(f_phis,req_i1)
109	CALL transfert_request(f_ps,req_i1)
110	CALL transfert_request(f_theta_rhodz,req_i1)
111	CALL transfert_request(f_u,req_e1)
112	CALL transfert_request(f_u,req_e1)
113
114
115	! CALL vorticity(f_u,f_out_z)
116
117	DO ind=1,ndomain
118	CALL swap_dimensions(ind)
119	CALL swap_geometry(ind)
120	CALL swap_caldyn(ind)
121
122	phis=f_phis(ind)
123	ps=f_ps(ind)
124	theta_rhodz=f_theta_rhodz(ind)
125	u=f_u(ind)
126	dps=f_dps(ind)
127	dtheta_rhodz=f_dtheta_rhodz(ind)
128	du=f_du(ind)
129
130	!$OMP PARALLEL DEFAULT(SHARED)
131	CALL compute_caldyn(phis, ps, theta_rhodz, u, dps, dtheta_rhodz, du)
132	!$OMP END PARALLEL
133	ENDDO
134
135	CALL transfert_request(f_out_u,req_e1)
136	CALL transfert_request(f_out_u,req_e1)
137
138	! CALL vorticity(f_u,f_out_z)
139	! CALL kinetic(f_du,f_out)
140
141	IF (mod(it,72)==0 ) THEN
142	CALL writefield("ps",f_ps)
143	CALL writefield("dps",f_dps)
144	! CALL writefield("theta_rhodz",f_theta_rhodz)
145	CALL writefield("dtheta_rhodz",f_dtheta_rhodz)
146	! CALL writefield("vort",f_out_z)
147	CALL writefield("theta",f_out)
148
149	! CALL writefield("out",f_out)
150	! DO ind=1,ndomain
151	! CALL writefield("Ki",f_out,ind)
152	! CALL writefield("vort",f_out_z,ind)
153	! ENDDO
154
155	ENDIF
156	! CALL check_mass_conservation(f_ps,f_dps)
157	it=it+1
158	END SUBROUTINE caldyn
159
160
161	SUBROUTINE compute_caldyn(phis, ps, theta_rhodz, u, dps, dtheta_rhodz, du)
162	USE domain_mod
163	USE dimensions
164	USE grid_param
165	USE geometry
166	USE metric
167	USE disvert_mod
168	IMPLICIT NONE
169	REAL(rstd),INTENT(IN) :: phis(iim*jjm)
170	REAL(rstd),INTENT(IN) :: u(iim3jjm,llm)
171	REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm)
172	REAL(rstd),INTENT(IN) :: ps(iim*jjm)
173	REAL(rstd),INTENT(OUT) :: du(iim3jjm,llm)
174	REAL(rstd),INTENT(OUT):: dtheta_rhodz(iim*jjm,llm)
175	REAL(rstd),INTENT(OUT):: dps(iim*jjm)
176
177	INTEGER :: i,j,ij,l
178	REAL(rstd) :: ww,uu
179	REAL(rstd) :: delta
180	REAL(rstd) :: etav,hv
181
182	! REAL(rstd) :: theta(iim*jjm,llm) ! potential temperature
183	! REAL(rstd) :: p(iim*jjm,llm+1) ! pression
184	! REAL(rstd) :: pk(iim*jjm,llm) ! Exner function
185	! REAL(rstd) :: pks(iim*jjm)
186	!! Intermediate variable to compute exner function
187	! REAL(rstd) :: alpha(iim*jjm,llm)
188	! REAL(rstd) :: beta(iim*jjm,llm)
189	!!
190	! REAL(rstd) :: phi(iim*jjm,llm) ! geopotential
191	! REAL(rstd) :: mass(iim*jjm,llm) ! mass
192	! REAL(rstd) :: rhodz(iim*jjm,llm) ! mass density
193	! REAL(rstd) :: Fe(3iimjjm,llm) ! mass flux
194	! REAL(rstd) :: Ftheta(3iimjjm,llm) ! theta flux
195	! REAL(rstd) :: convm(iim*jjm,llm) ! mass flux convergence
196	! REAL(rstd) :: w(iim*jjm,llm) ! vertical velocity
197	! REAL(rstd) :: qv(2iimjjm,llm) ! potential velocity
198	! REAL(rstd) :: berni(iim*jjm,llm) ! bernouilli term
199
200	REAL(rstd),ALLOCATABLE,SAVE :: theta(:,:) ! potential temperature
201	REAL(rstd),ALLOCATABLE,SAVE :: p(:,:) ! pression
202	REAL(rstd),ALLOCATABLE,SAVE :: pk(:,:) ! Exner function
203	REAL(rstd),ALLOCATABLE,SAVE :: pks(:)
204	! Intermediate variable to compute exner function
205	REAL(rstd),ALLOCATABLE,SAVE :: alpha(:,:)
206	REAL(rstd),ALLOCATABLE,SAVE :: beta(:,:)
207	!
208	REAL(rstd),ALLOCATABLE,SAVE :: phi(:,:) ! geopotential
209	REAL(rstd),ALLOCATABLE,SAVE :: mass(:,:) ! mass
210	REAL(rstd),ALLOCATABLE,SAVE :: rhodz(:,:) ! mass density
211	REAL(rstd),ALLOCATABLE,SAVE :: Fe(:,:) ! mass flux
212	REAL(rstd),ALLOCATABLE,SAVE :: Ftheta(:,:) ! theta flux
213	REAL(rstd),ALLOCATABLE,SAVE :: convm(:,:) ! mass flux convergence
214	REAL(rstd),ALLOCATABLE,SAVE :: w(:,:) ! vertical velocity
215	REAL(rstd),ALLOCATABLE,SAVE :: qv(:,:) ! potential velocity
216	REAL(rstd),ALLOCATABLE,SAVE :: berni(:,:) ! bernouilli term
217
218	LOGICAL,SAVE :: first=.TRUE.
219	!$OMP THREADPRIVATE(first)
220
221	!$OMP BARRIER
222	!$OMP MASTER
223	IF (first) THEN
224	ALLOCATE(theta(iim*jjm,llm)) ! potential temperature
225	ALLOCATE(p(iim*jjm,llm+1)) ! pression
226	ALLOCATE(pk(iim*jjm,llm)) ! Exner function
227	ALLOCATE(pks(iim*jjm))
228	ALLOCATE(alpha(iim*jjm,llm))
229	ALLOCATE(beta(iim*jjm,llm))
230	ALLOCATE(phi(iim*jjm,llm)) ! geopotential
231	ALLOCATE(mass(iim*jjm,llm)) ! mass
232	ALLOCATE(rhodz(iim*jjm,llm)) ! mass density
233	ALLOCATE(Fe(3iimjjm,llm)) ! mass flux
234	ALLOCATE(Ftheta(3iimjjm,llm)) ! theta flux
235	ALLOCATE(convm(iim*jjm,llm)) ! mass flux convergence
236	ALLOCATE(w(iim*jjm,llm)) ! vertical velocity
237	ALLOCATE(qv(2iimjjm,llm)) ! potential velocity
238	ALLOCATE(berni(iim*jjm,llm)) ! bernouilli term
239	first=.FALSE.
240	ENDIF
241	!$OMP END MASTER
242	!$OMP BARRIER
243	! du(:,:)=0
244	! theta=1e10
245	! p=1e10
246	! pk=1e10
247	! pks=1e10
248	! alpha=1e10
249	! beta=1e10
250	! phi=1e10
251	! mass=1e10
252	! rhodz=1e10
253	! Fe=1e10
254	! Ftheta=1e10
255	! convm=1e10
256	! w=1e10
257	! qv=1e10
258	! berni=1e10
259
260	!!! Compute pression
261
262	DO l = 1, llm+1
263	!$OMP DO
264	DO j=jj_begin-1,jj_end+1
265	DO i=ii_begin-1,ii_end+1
266	ij=(j-1)*iim+i
267	p(ij,l) = ap(l) + bp(l) * ps(ij)
268	ENDDO
269	ENDDO
270	ENDDO
271
272
273
274
275	!!! Compute Exnher function
276
277	!! Compute Alpha and Beta
278
279	! for llm layer
280	!$OMP DO
281	DO j=jj_begin-1,jj_end+1
282	DO i=ii_begin-1,ii_end+1
283	ij=(j-1)*iim+i
284	alpha(ij,llm) = 0.
285	beta (ij,llm) = 1./ (1+ 2*kappa)
286	ENDDO
287	ENDDO
288
289	! for other layer
290	DO l = llm-1 , 2 , -1
291	!$OMP DO
292	DO j=jj_begin-1,jj_end+1
293	DO i=ii_begin-1,ii_end+1
294	ij=(j-1)*iim+i
295	delta = p(ij,l)* (1+2kappa) + p(ij,l+1) ( beta(ij,l+1)- (1+2*kappa) )
296	alpha(ij,l) = - p(ij,l+1) / delta * alpha(ij,l+1)
297	beta (ij,l) = p(ij,l ) / delta
298	ENDDO
299	ENDDO
300	ENDDO
301
302	!! Compute pk
303
304	! for first layer
305	!$OMP DO
306	DO j=jj_begin-1,jj_end+1
307	DO i=ii_begin-1,ii_end+1
308	ij=(j-1)*iim+i
309	pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
310	pk(ij,1) = ( p(ij,1)pks(ij) - 0.5alpha(ij,2)*p(ij,2) ) / &
311	( p(ij,1)* (1.+kappa) + 0.5( beta(ij,2)-(1.+2kappa) )* p(ij,2) )
312	ENDDO
313	ENDDO
314
315	! for other layers
316
317	DO l = 2, llm
318	!$OMP DO
319	DO j=jj_begin-1,jj_end+1
320	DO i=ii_begin-1,ii_end+1
321	ij=(j-1)*iim+i
322	pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
323	ENDDO
324	ENDDO
325	ENDDO
326
327
328	!!! Compute mass
329	DO l = 1, llm
330	!$OMP DO
331	DO j=jj_begin-1,jj_end+1
332	DO i=ii_begin-1,ii_end+1
333	ij=(j-1)*iim+i
334	mass(ij,l) = ( p(ij,l) - p(ij,l+1) ) * Ai(ij)/g
335	rhodz(ij,l) = mass(ij,l) / Ai(ij)
336	ENDDO
337	ENDDO
338	ENDDO
339
340	!! compute theta
341	DO l = 1, llm
342	!$OMP DO
343	DO j=jj_begin-1,jj_end+1
344	DO i=ii_begin-1,ii_end+1
345	ij=(j-1)*iim+i
346	theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
347	ENDDO
348	ENDDO
349	ENDDO
350
351
352	!!! Compute geopotential
353
354	! for first layer
355	!$OMP DO
356	DO j=jj_begin-1,jj_end+1
357	DO i=ii_begin-1,ii_end+1
358	ij=(j-1)*iim+i
359	phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
360	ENDDO
361	ENDDO
362
363	! for other layers
364	DO l = 2, llm
365	!$OMP DO
366	DO j=jj_begin-1,jj_end+1
367	DO i=ii_begin-1,ii_end+1
368	ij=(j-1)*iim+i
369	phi(ij,l) = phi(ij,l-1) + 0.5 * ( theta(ij,l) + theta(ij,l-1) ) &
370	* ( pk(ij,l-1) - pk(ij,l) )
371	ENDDO
372	ENDDO
373	ENDDO
374
375
376	!!! Compute mass flux
377	!! question à thomas : meilleure pondération de la masse sur les liens ?
378
379	DO l = 1, llm
380	!$OMP DO
381	DO j=jj_begin-1,jj_end+1
382	DO i=ii_begin-1,ii_end+1
383	ij=(j-1)*iim+i
384	Fe(ij+u_right,l)=0.5(rhodz(ij,l)+rhodz(ij+t_right,l))u(ij+u_right,l)*le(ij+u_right)
385	Fe(ij+u_lup,l)=0.5(rhodz(ij,l)+rhodz(ij+t_lup,l))u(ij+u_lup,l)*le(ij+u_lup)
386	Fe(ij+u_ldown,l)=0.5(rhodz(ij,l)+rhodz(ij+t_ldown,l))u(ij+u_ldown,l)*le(ij+u_ldown)
387	ENDDO
388	ENDDO
389	ENDDO
390
391	!!! fisrt composante dtheta
392
393	! Flux on the edge
394	DO l = 1, llm
395	!$OMP DO
396	DO j=jj_begin-1,jj_end+1
397	DO i=ii_begin-1,ii_end+1
398	ij=(j-1)*iim+i
399	Ftheta(ij+u_right,l)=0.5(theta(ij,l)+theta(ij+t_right,l))Fe(ij+u_right,l)
400	Ftheta(ij+u_lup,l)=0.5(theta(ij,l)+theta(ij+t_lup,l))Fe(ij+u_lup,l)
401	Ftheta(ij+u_ldown,l)=0.5(theta(ij,l)+theta(ij+t_ldown,l))Fe(ij+u_ldown,l)
402	ENDDO
403	ENDDO
404	ENDDO
405
406
407	! compute divergence
408	DO l = 1, llm
409	!$OMP DO
410	DO j=jj_begin,jj_end
411	DO i=ii_begin,ii_end
412	ij=(j-1)*iim+i
413	! signe ? attention d (rho theta dz)
414	dtheta_rhodz(ij,l)=-1./Ai(ij)(ne(ij,right)Ftheta(ij+u_right,l) + &
415	ne(ij,rup)*Ftheta(ij+u_rup,l) + &
416	ne(ij,lup)*Ftheta(ij+u_lup,l) + &
417	ne(ij,left)*Ftheta(ij+u_left,l) + &
418	ne(ij,ldown)*Ftheta(ij+u_ldown,l) + &
419	ne(ij,rdown)*Ftheta(ij+u_rdown,l))
420	ENDDO
421	ENDDO
422	ENDDO
423
424
425
426	!!! mass flux convergence computation
427
428	! horizontal convergence
429	DO l = 1, llm
430	!$OMP DO
431	DO j=jj_begin,jj_end
432	DO i=ii_begin,ii_end
433	ij=(j-1)*iim+i
434	!signe ?
435	convm(ij,l)= 1./Ai(ij)(ne(ij,right)Fe(ij+u_right,l) + &
436	ne(ij,rup)*Fe(ij+u_rup,l) + &
437	ne(ij,lup)*Fe(ij+u_lup,l) + &
438	ne(ij,left)*Fe(ij+u_left,l) + &
439	ne(ij,ldown)*Fe(ij+u_ldown,l) + &
440	ne(ij,rdown)*Fe(ij+u_rdown,l))
441	ENDDO
442	ENDDO
443	ENDDO
444
445
446	! vertical integration from up to down
447	DO l = llm-1, 1, -1
448	!$OMP DO
449	DO j=jj_begin,jj_end
450	DO i=ii_begin,ii_end
451	ij=(j-1)*iim+i
452	convm(ij,l) = convm(ij,l) + convm(ij,l+1)
453	ENDDO
454	ENDDO
455	ENDDO
456
457
458	DO l = 1, llm
459	!$OMP DO
460	DO j=jj_begin,jj_end
461	DO i=ii_begin,ii_end
462	ij=(j-1)*iim+i
463	out(ij,l)=theta(ij,l)-288
464	ENDDO
465	ENDDO
466	ENDDO
467
468
469	!!! Compute dps
470	!$OMP DO
471	DO j=jj_begin,jj_end
472	DO i=ii_begin,ii_end
473	ij=(j-1)*iim+i
474	dps(ij)=-convm(ij,1) * g
475	ENDDO
476	ENDDO
477
478
479
480	!!! Compute vertical velocity
481	DO l = 1,llm-1
482	!$OMP DO
483	DO j=jj_begin,jj_end
484	DO i=ii_begin,ii_end
485	ij=(j-1)*iim+i
486	w( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 )
487	ENDDO
488	ENDDO
489	ENDDO
490
491	!$OMP DO
492	DO j=jj_begin,jj_end
493	DO i=ii_begin,ii_end
494	ij=(j-1)*iim+i
495	w(ij,1) = 0.
496	ENDDO
497	ENDDO
498
499
500	!!! Compute potential vorticity
501	DO l = 1,llm
502	!$OMP DO
503	DO j=jj_begin-1,jj_end+1
504	DO i=ii_begin-1,ii_end+1
505	ij=(j-1)*iim+i
506
507	etav= 1./Av(ij+z_up)( ne(ij,rup) u(ij+u_rup,l) * de(ij+u_rup) &
508	+ ne(ij+t_rup,left) * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left) &
509	- ne(ij,lup) * u(ij+u_lup,l) * de(ij+u_lup) )
510
511	hv = Riv2(ij,vup) * rhodz(ij,l) &
512	+ Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l) &
513	+ Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
514
515	qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv
516
517	etav = 1./Av(ij+z_down)( ne(ij,ldown) u(ij+u_ldown,l) * de(ij+u_ldown) &
518	+ ne(ij+t_ldown,right) * u(ij+t_ldown+u_right,l) * de(ij+t_ldown+u_right) &
519	- ne(ij,rdown) * u(ij+u_rdown,l) * de(ij+u_rdown) )
520
521	hv = Riv2(ij,vdown) * rhodz(ij,l) &
522	+ Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l) &
523	+ Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
524
525	qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv
526
527	ENDDO
528	ENDDO
529	ENDDO
530
531	!!! Compute potential vorticity contribution to du
532	DO l=1,llm
533	!$OMP DO
534	DO j=jj_begin,jj_end
535	DO i=ii_begin,ii_end
536	ij=(j-1)*iim+i
537
538	du(ij+u_right,l) = 0.5(qv(ij+z_rdown,l)+qv(ij+z_rup,l))/de(ij+u_right) &
539	( wee(ij+u_right,1,1)*Fe(ij+u_rup,l)+ &
540	wee(ij+u_right,2,1)*Fe(ij+u_lup,l)+ &
541	wee(ij+u_right,3,1)*Fe(ij+u_left,l)+ &
542	wee(ij+u_right,4,1)*Fe(ij+u_ldown,l)+ &
543	wee(ij+u_right,5,1)*Fe(ij+u_rdown,l)+ &
544	wee(ij+u_right,1,2)*Fe(ij+t_right+u_ldown,l)+ &
545	wee(ij+u_right,2,2)*Fe(ij+t_right+u_rdown,l)+ &
546	wee(ij+u_right,3,2)*Fe(ij+t_right+u_right,l)+ &
547	wee(ij+u_right,4,2)*Fe(ij+t_right+u_rup,l)+ &
548	wee(ij+u_right,5,2)*Fe(ij+t_right+u_lup,l) )
549
550
551	du(ij+u_lup,l) = 0.5(qv(ij+z_up,l)+qv(ij+z_lup,l))/de(ij+u_lup) &
552	( wee(ij+u_lup,1,1)*Fe(ij+u_left,l)+ &
553	wee(ij+u_lup,2,1)*Fe(ij+u_ldown,l)+ &
554	wee(ij+u_lup,3,1)*Fe(ij+u_rdown,l)+ &
555	wee(ij+u_lup,4,1)*Fe(ij+u_right,l)+ &
556	wee(ij+u_lup,5,1)*Fe(ij+u_rup,l)+ &
557	wee(ij+u_lup,1,2)*Fe(ij+t_lup+u_right,l)+ &
558	wee(ij+u_lup,2,2)*Fe(ij+t_lup+u_rup,l)+ &
559	wee(ij+u_lup,3,2)*Fe(ij+t_lup+u_lup,l)+ &
560	wee(ij+u_lup,4,2)*Fe(ij+t_lup+u_left,l)+ &
561	wee(ij+u_lup,5,2)*Fe(ij+t_lup+u_ldown,l) )
562
563
564	du(ij+u_ldown,l) = 0.5(qv(ij+z_ldown,l)+qv(ij+z_down,l))/de(ij+u_ldown) &
565	( wee(ij+u_ldown,1,1)*Fe(ij+u_rdown,l)+ &
566	wee(ij+u_ldown,2,1)*Fe(ij+u_right,l)+ &
567	wee(ij+u_ldown,3,1)*Fe(ij+u_rup,l)+ &
568	wee(ij+u_ldown,4,1)*Fe(ij+u_lup,l)+ &
569	wee(ij+u_ldown,5,1)*Fe(ij+u_left,l)+ &
570	wee(ij+u_ldown,1,2)*Fe(ij+t_ldown+u_lup,l)+ &
571	wee(ij+u_ldown,2,2)*Fe(ij+t_ldown+u_left,l)+ &
572	wee(ij+u_ldown,3,2)*Fe(ij+t_ldown+u_ldown,l)+ &
573	wee(ij+u_ldown,4,2)*Fe(ij+t_ldown+u_rdown,l)+ &
574	wee(ij+u_ldown,5,2)*Fe(ij+t_ldown+u_right,l) )
575
576
577	ENDDO
578	ENDDO
579	ENDDO
580
581
582	!!! Compute bernouilli term = Kinetic Energy + geopotential
583	DO l=1,llm
584	!$OMP DO
585	DO j=jj_begin,jj_end
586	DO i=ii_begin,ii_end
587	ij=(j-1)*iim+i
588
589	berni(ij,l) = phi(ij,l) &
590	+ 1/(4Ai(ij))(le(ij+u_right)de(ij+u_right)u(ij+u_right,l)**2 + &
591	le(ij+u_rup)de(ij+u_rup)u(ij+u_rup,l)**2 + &
592	le(ij+u_lup)de(ij+u_lup)u(ij+u_lup,l)**2 + &
593	le(ij+u_left)de(ij+u_left)u(ij+u_left,l)**2 + &
594	le(ij+u_ldown)de(ij+u_ldown)u(ij+u_ldown,l)**2 + &
595	le(ij+u_rdown)de(ij+u_rdown)u(ij+u_rdown,l)**2 )
596
597	ENDDO
598	ENDDO
599	ENDDO
600
601
602	!!! second contribution to du
603	DO l=1,llm
604	!$OMP DO
605	DO j=jj_begin,jj_end
606	DO i=ii_begin,ii_end
607	ij=(j-1)*iim+i
608
609	du(ij+u_right,l)= du(ij+u_right,l)+ 1/de(ij+u_right) * ( &
610	0.5*(theta(ij,l)+theta(ij+t_right,l)) &
611	( ne(ij,right)pk(ij,l)+ne(ij+t_right,left)*pk(ij+t_right,l)) &
612	+ ne(ij,right)berni(ij,l)+ne(ij+t_right,left)berni(ij+t_right,l) )
613
614	du(ij+u_lup,l)= du(ij+u_lup,l)+ 1/de(ij+u_lup) * ( &
615	0.5*(theta(ij,l)+theta(ij+t_lup,l)) &
616	( ne(ij,lup)pk(ij,l)+ne(ij+t_lup,rdown)*pk(ij+t_lup,l)) &
617	+ ne(ij,lup)berni(ij,l)+ne(ij+t_lup,rdown)berni(ij+t_lup,l) )
618
619	du(ij+u_ldown,l)= du(ij+u_ldown,l)+ 1/de(ij+u_ldown) * ( &
620	0.5*(theta(ij,l)+theta(ij+t_ldown,l)) &
621	( ne(ij,ldown)pk(ij,l)+ne(ij+t_ldown,rup)*pk(ij+t_ldown,l)) &
622	+ ne(ij,ldown)berni(ij,l)+ne(ij+t_ldown,rup)berni(ij+t_ldown,l) )
623	ENDDO
624	ENDDO
625	ENDDO
626
627	!!! second contribution to du
628	DO l=1,llm
629	!$OMP DO
630	DO j=jj_begin,jj_end
631	DO i=ii_begin,ii_end
632	ij=(j-1)*iim+i
633
634	out_u(ij+u_right,l)= 1/de(ij+u_right) * ( &
635	0.5*(theta(ij,l)+theta(ij+t_right,l)) &
636	( ne(ij,right)pk(ij,l)+ne(ij+t_right,left)*pk(ij+t_right,l)) &
637	+ ne(ij,right)berni(ij,l)+ne(ij+t_right,left)berni(ij+t_right,l) )
638
639	out_u(ij+u_lup,l)= 1/de(ij+u_lup) * ( &
640	0.5*(theta(ij,l)+theta(ij+t_lup,l)) &
641	( ne(ij,lup)pk(ij,l)+ne(ij+t_lup,rdown)*pk(ij+t_lup,l)) &
642	+ ne(ij,lup)berni(ij,l)+ne(ij+t_lup,rdown)berni(ij+t_lup,l) )
643
644	out_u(ij+u_ldown,l)= 1/de(ij+u_ldown) * ( &
645	0.5*(theta(ij,l)+theta(ij+t_ldown,l)) &
646	( ne(ij,ldown)pk(ij,l)+ne(ij+t_ldown,rup)*pk(ij+t_ldown,l)) &
647	+ ne(ij,ldown)berni(ij,l)+ne(ij+t_ldown,rup)berni(ij+t_ldown,l) )
648	ENDDO
649	ENDDO
650	ENDDO
651	!!! contribution due to vertical advection
652
653	! Contribution to dtheta
654	DO l=1,llm-1
655	!$OMP DO
656	DO j=jj_begin,jj_end
657	DO i=ii_begin,ii_end
658	ij=(j-1)*iim+i
659	ww = 0.5 * w(ij,l+1) * (theta(ij,l) + theta(ij,l+1) )
660	dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l ) - ww
661	dtheta_rhodz(ij,l+1) = dtheta_rhodz(ij,l+1) + ww
662	ENDDO
663	ENDDO
664	ENDDO
665
666
667	! Contribution to du
668	DO l=1,llm-1
669	!$OMP DO
670	DO j=jj_begin,jj_end
671	DO i=ii_begin,ii_end
672	ij=(j-1)*iim+i
673	ww = 0.5 * ( w(ij,l+1) + w(ij+t_right,l+1))
674	uu = u(ij+u_right,l+1) - u(ij+u_right,l)
675	du(ij+u_right, l ) = du(ij+u_right,l) - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_right,l)))
676	du(ij+u_right, l+1 ) = du(ij+u_right,l+1) - 0.5 * ww * uu / (0.5*(rhodz(ij,l+1)+rhodz(ij+t_right,l+1)))
677
678	ww = 0.5 * ( w(ij,l+1) + w(ij+t_lup,l+1))
679	uu = u(ij+u_lup,l+1) - u(ij+u_lup,l)
680	du(ij+u_lup, l ) = du(ij+u_lup,l) - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l)))
681	du(ij+u_lup, l+1 ) = du(ij+u_lup,l+1) - 0.5 * ww * uu / (0.5*(rhodz(ij,l+1)+rhodz(ij+t_lup,l+1)))
682
683	ww = 0.5 * ( w(ij,l+1) + w(ij+t_ldown,l+1))
684	uu = u(ij+u_ldown,l+1) - u(ij+u_ldown,l)
685	du(ij+u_ldown, l ) = du(ij+u_ldown,l) - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l)))
686	du(ij+u_ldown, l+1 ) = du(ij+u_ldown,l+1) - 0.5 * ww * uu / (0.5*(rhodz(ij,l+1)+rhodz(ij+t_ldown,l+1)))
687
688	ENDDO
689	ENDDO
690	ENDDO
691
692	!!$OMP BARRIER
693	!!$OMP MASTER
694	! DEALLOCATE(theta) ! potential temperature
695	! DEALLOCATE(p) ! pression
696	! DEALLOCATE(pk) ! Exner function
697	! DEALLOCATE(pks)
698	! DEALLOCATE(alpha)
699	! DEALLOCATE(beta)
700	! DEALLOCATE(phi) ! geopotential
701	! DEALLOCATE(mass) ! mass
702	! DEALLOCATE(rhodz) ! mass density
703	! DEALLOCATE(Fe) ! mass flux
704	! DEALLOCATE(Ftheta) ! theta flux
705	! DEALLOCATE(convm) ! mass flux convergence
706	! DEALLOCATE(w) ! vertical velocity
707	! DEALLOCATE(qv) ! potential velocity
708	! DEALLOCATE(berni) ! bernouilli term
709	!!$OMP END MASTER
710	!!$OMP BARRIER
711	END SUBROUTINE compute_caldyn
712
713
714
715	END MODULE caldyn_gcm_mod

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

Download in other formats: