source: codes/icosagcm/trunk/src/caldyn_adv.f90 @ 111

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

Correction for dcmip moist physics
Temperature is ouput instead of virtual temperature.

YM
File size: 7.9 KB
Line 
1MODULE caldyn_adv_mod
2  USE icosa
3
4  TYPE(t_field),POINTER :: f_out(:)
5  REAL(rstd),POINTER :: out(:,:)
6  TYPE(t_field),POINTER :: f_out_u(:)
7  REAL(rstd),POINTER :: out_u(:,:)
8  TYPE(t_field),POINTER :: f_out_z(:)
9  REAL(rstd),POINTER :: out_z(:,:)
10 
11 
12CONTAINS
13 
14  SUBROUTINE init_caldyn
15  USE icosa
16  IMPLICIT NONE
17
18    CALL allocate_caldyn
19 
20  END SUBROUTINE init_caldyn
21 
22  SUBROUTINE allocate_caldyn
23  USE icosa
24  IMPLICIT NONE
25
26    CALL allocate_field(f_out,field_t,type_real,llm) 
27    CALL allocate_field(f_out_u,field_u,type_real,llm) 
28    CALL allocate_field(f_out_z,field_z,type_real,llm) 
29   
30  END SUBROUTINE allocate_caldyn
31 
32  SUBROUTINE swap_caldyn(ind)
33  IMPLICIT NONE
34    INTEGER,INTENT(IN) :: ind
35    out=f_out(ind) 
36    out_u=f_out_u(ind) 
37    out_z=f_out_z(ind) 
38     
39  END SUBROUTINE swap_caldyn
40 
41  SUBROUTINE check_mass_conservation(f_ps,f_dps)
42  USE icosa
43  IMPLICIT NONE
44    TYPE(t_field),POINTER :: f_ps(:)
45    TYPE(t_field),POINTER :: f_dps(:)
46    REAL(rstd),POINTER :: ps(:)
47    REAL(rstd),POINTER :: dps(:)
48    REAL(rstd) :: mass_tot,dmass_tot
49    INTEGER :: ind,i,j,ij
50   
51    mass_tot=0
52    dmass_tot=0
53   
54    CALL transfert_request(f_dps,req_i1)
55    CALL transfert_request(f_ps,req_i1)
56
57    DO ind=1,ndomain
58      CALL swap_dimensions(ind)
59      CALL swap_geometry(ind)
60
61      ps=f_ps(ind)
62      dps=f_dps(ind)
63
64      DO j=jj_begin,jj_end
65        DO i=ii_begin,ii_end
66          ij=(j-1)*iim+i
67          IF (domain(ind)%own(i,j)) THEN
68            mass_tot=mass_tot+ps(ij)*Ai(ij)/g
69            dmass_tot=dmass_tot+dps(ij)*Ai(ij)/g
70          ENDIF
71        ENDDO
72      ENDDO
73   
74    ENDDO
75    PRINT*, "mass_tot ", mass_tot,"      dmass_tot ",dmass_tot       
76
77  END SUBROUTINE check_mass_conservation
78 
79 
80 
81  SUBROUTINE caldyn(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du)
82  USE icosa
83  USE vorticity_mod
84  USE kinetic_mod
85  USE theta2theta_rhodz_mod
86  IMPLICIT NONE
87  INTEGER,INTENT(IN)    :: it
88  TYPE(t_field),POINTER :: f_phis(:)
89  TYPE(t_field),POINTER :: f_ps(:)
90  TYPE(t_field),POINTER :: f_theta_rhodz(:)
91  TYPE(t_field),POINTER :: f_u(:)
92  TYPE(t_field),POINTER :: f_q(:)
93  TYPE(t_field),POINTER :: f_dps(:)
94  TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
95  TYPE(t_field),POINTER :: f_du(:)
96
97  REAL(rstd),POINTER :: phis(:)
98  REAL(rstd),POINTER :: ps(:)
99  REAL(rstd),POINTER :: theta_rhodz(:,:)
100  REAL(rstd),POINTER :: u(:,:)
101  REAL(rstd),POINTER :: dps(:)
102  REAL(rstd),POINTER :: dtheta_rhodz(:,:)
103  REAL(rstd),POINTER :: du(:,:)
104  INTEGER :: ind
105
106 
107    CALL transfert_request(f_phis,req_i1) 
108    CALL transfert_request(f_ps,req_i1) 
109    CALL transfert_request(f_theta_rhodz,req_i1) 
110    CALL transfert_request(f_u,req_e1)
111    CALL transfert_request(f_u,req_e1) 
112 
113    DO ind=1,ndomain
114      CALL swap_dimensions(ind)
115      CALL swap_geometry(ind)
116      CALL swap_caldyn(ind)
117     
118      phis=f_phis(ind)
119      ps=f_ps(ind)
120      theta_rhodz=f_theta_rhodz(ind)
121      u=f_u(ind)
122      dps=f_dps(ind)
123      dtheta_rhodz=f_dtheta_rhodz(ind)
124      du=f_du(ind)
125     
126!$OMP PARALLEL DEFAULT(SHARED)
127      CALL compute_caldyn(phis, ps, theta_rhodz, u, dps, dtheta_rhodz, du)
128!$OMP END PARALLEL
129    ENDDO
130
131    CALL transfert_request(f_out_u,req_e1)
132    CALL transfert_request(f_out_u,req_e1) 
133
134
135    IF (mod(it,itau_out)==0 ) THEN
136      CALL writefield("ps",f_ps)
137    ENDIF
138!    CALL check_mass_conservation(f_ps,f_dps)
139
140  END SUBROUTINE caldyn
141 
142 
143  SUBROUTINE compute_caldyn(phis, ps, theta_rhodz, u, dps, dtheta_rhodz, du)
144  USE icosa
145  USE disvert_mod
146  IMPLICIT NONE
147    REAL(rstd),INTENT(IN) :: phis(iim*jjm)
148    REAL(rstd),INTENT(IN) :: u(iim*3*jjm,llm)
149    REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm)
150    REAL(rstd),INTENT(IN) :: ps(iim*jjm)
151    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm)
152    REAL(rstd),INTENT(OUT):: dtheta_rhodz(iim*jjm,llm)
153    REAL(rstd),INTENT(OUT):: dps(iim*jjm)
154   
155    INTEGER :: i,j,ij,l
156    REAL(rstd) :: ww,uu 
157    REAL(rstd) :: delta
158    REAL(rstd) :: etav,hv   
159    REAL(rstd),ALLOCATABLE,SAVE :: theta(:,:)  ! potential temperature
160    REAL(rstd),ALLOCATABLE,SAVE :: p(:,:)  ! pression
161    REAL(rstd),ALLOCATABLE,SAVE :: pk(:,:)   ! Exner function
162    REAL(rstd),ALLOCATABLE,SAVE :: pks(:)
163 ! Intermediate variable to compute exner function
164    REAL(rstd),ALLOCATABLE,SAVE :: alpha(:,:)
165    REAL(rstd),ALLOCATABLE,SAVE :: beta(:,:)
166 !   
167    REAL(rstd),ALLOCATABLE,SAVE :: phi(:,:)    ! geopotential
168    REAL(rstd),ALLOCATABLE,SAVE :: mass(:,:)   ! mass   
169    REAL(rstd),ALLOCATABLE,SAVE :: rhodz(:,:)   ! mass density   
170    REAL(rstd),ALLOCATABLE,SAVE :: Fe(:,:)   ! mass flux   
171    REAL(rstd),ALLOCATABLE,SAVE :: Ftheta(:,:) ! theta flux   
172    REAL(rstd),ALLOCATABLE,SAVE :: convm(:,:)    ! mass flux convergence
173    REAL(rstd),ALLOCATABLE,SAVE :: w(:,:)       ! vertical velocity     
174    REAL(rstd),ALLOCATABLE,SAVE :: qv(:,:)       ! potential velocity 
175    REAL(rstd),ALLOCATABLE,SAVE :: berni(:,:)   ! bernouilli term
176   
177    LOGICAL,SAVE :: first=.TRUE.
178!$OMP THREADPRIVATE(first)
179       
180!$OMP BARRIER     
181!$OMP MASTER 
182    IF (first) THEN
183      ALLOCATE(theta(iim*jjm,llm))  ! potential temperature
184      ALLOCATE(p(iim*jjm,llm+1))  ! pression
185      ALLOCATE(pk(iim*jjm,llm))   ! Exner function
186      ALLOCATE(pks(iim*jjm))
187      ALLOCATE(alpha(iim*jjm,llm))
188      ALLOCATE(beta(iim*jjm,llm))
189      ALLOCATE(phi(iim*jjm,llm))    ! geopotential
190      ALLOCATE(mass(iim*jjm,llm))   ! mass   
191      ALLOCATE(rhodz(iim*jjm,llm))   ! mass density   
192      ALLOCATE(Fe(3*iim*jjm,llm))   ! mass flux   
193      ALLOCATE(Ftheta(3*iim*jjm,llm)) ! theta flux   
194      ALLOCATE(convm(iim*jjm,llm))    ! mass flux convergence
195      ALLOCATE(w(iim*jjm,llm))       ! vertical velocity     
196      ALLOCATE(qv(2*iim*jjm,llm))       ! potential velocity 
197      ALLOCATE(berni(iim*jjm,llm))   ! bernouilli term
198      first=.FALSE.
199    ENDIF
200!$OMP END MASTER
201!$OMP BARRIER       
202   
203 !!! Compute pression
204   dtheta_rhodz(:,:)=0.
205   du(:,:)=0.
206 
207    DO    l    = 1, llm+1
208!$OMP DO
209      DO j=jj_begin-1,jj_end+1
210        DO i=ii_begin-1,ii_end+1
211          ij=(j-1)*iim+i
212          p(ij,l) = ap(l) + bp(l) * ps(ij)
213        ENDDO
214      ENDDO
215    ENDDO
216     
217 
218!!! Compute mass
219   DO l = 1, llm
220!$OMP DO
221     DO j=jj_begin-1,jj_end+1
222       DO i=ii_begin-1,ii_end+1
223         ij=(j-1)*iim+i
224         mass(ij,l) = ( p(ij,l) - p(ij,l+1) ) * Ai(ij)/g
225         rhodz(ij,l) = mass(ij,l) / Ai(ij)
226       ENDDO
227     ENDDO
228   ENDDO
229
230   
231!!!  Compute mass flux
232!! question à thomas : meilleure pondération de la masse sur les liens ?
233
234  DO l = 1, llm
235!$OMP DO
236    DO j=jj_begin-1,jj_end+1
237      DO i=ii_begin-1,ii_end+1
238        ij=(j-1)*iim+i
239        Fe(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)
240        Fe(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)
241        Fe(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown)
242      ENDDO
243    ENDDO
244  ENDDO
245
246
247!!! mass flux convergence computation 
248
249 ! horizontal convergence
250  DO l = 1, llm
251!$OMP DO
252    DO j=jj_begin,jj_end
253      DO i=ii_begin,ii_end
254        ij=(j-1)*iim+i
255!signe ?
256        convm(ij,l)= 1./Ai(ij)*(ne(ij,right)*Fe(ij+u_right,l)   +  &
257                                ne(ij,rup)*Fe(ij+u_rup,l)       +  & 
258                                ne(ij,lup)*Fe(ij+u_lup,l)       +  & 
259                                ne(ij,left)*Fe(ij+u_left,l)     +  & 
260                                ne(ij,ldown)*Fe(ij+u_ldown,l)   +  & 
261                                ne(ij,rdown)*Fe(ij+u_rdown,l))   
262       ENDDO
263     ENDDO
264   ENDDO
265 
266   
267 ! vertical integration from up to down
268  DO  l = llm-1, 1, -1
269!$OMP DO
270    DO j=jj_begin,jj_end
271      DO i=ii_begin,ii_end
272        ij=(j-1)*iim+i
273        convm(ij,l) = convm(ij,l) + convm(ij,l+1)
274      ENDDO
275    ENDDO
276  ENDDO
277 
278
279!!! Compute dps
280!$OMP DO
281  DO j=jj_begin,jj_end
282    DO i=ii_begin,ii_end
283      ij=(j-1)*iim+i
284      dps(ij)=-convm(ij,1) * g 
285    ENDDO
286  ENDDO
287                           
288  END SUBROUTINE compute_caldyn
289 
290 
291 
292END MODULE caldyn_adv_mod
Note: See TracBrowser for help on using the repository browser.