source: codes/icosagcm/trunk/src/caldyn_adv.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: 7.8 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_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_dps(:)
93  TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
94  TYPE(t_field),POINTER :: f_du(:)
95
96  REAL(rstd),POINTER :: phis(:)
97  REAL(rstd),POINTER :: ps(:)
98  REAL(rstd),POINTER :: theta_rhodz(:,:)
99  REAL(rstd),POINTER :: u(:,:)
100  REAL(rstd),POINTER :: dps(:)
101  REAL(rstd),POINTER :: dtheta_rhodz(:,:)
102  REAL(rstd),POINTER :: du(:,:)
103  INTEGER :: ind
104
105 
106    CALL transfert_request(f_phis,req_i1) 
107    CALL transfert_request(f_ps,req_i1) 
108    CALL transfert_request(f_theta_rhodz,req_i1) 
109    CALL transfert_request(f_u,req_e1)
110    CALL transfert_request(f_u,req_e1) 
111 
112    DO ind=1,ndomain
113      CALL swap_dimensions(ind)
114      CALL swap_geometry(ind)
115      CALL swap_caldyn(ind)
116     
117      phis=f_phis(ind)
118      ps=f_ps(ind)
119      theta_rhodz=f_theta_rhodz(ind)
120      u=f_u(ind)
121      dps=f_dps(ind)
122      dtheta_rhodz=f_dtheta_rhodz(ind)
123      du=f_du(ind)
124     
125!$OMP PARALLEL DEFAULT(SHARED)
126      CALL compute_caldyn(phis, ps, theta_rhodz, u, dps, dtheta_rhodz, du)
127!$OMP END PARALLEL
128    ENDDO
129
130    CALL transfert_request(f_out_u,req_e1)
131    CALL transfert_request(f_out_u,req_e1) 
132
133
134    IF (mod(it,itau_out)==0 ) THEN
135      CALL writefield("ps",f_ps)
136    ENDIF
137!    CALL check_mass_conservation(f_ps,f_dps)
138
139  END SUBROUTINE caldyn
140 
141 
142  SUBROUTINE compute_caldyn(phis, ps, theta_rhodz, u, dps, dtheta_rhodz, du)
143  USE icosa
144  USE disvert_mod
145  IMPLICIT NONE
146    REAL(rstd),INTENT(IN) :: phis(iim*jjm)
147    REAL(rstd),INTENT(IN) :: u(iim*3*jjm,llm)
148    REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm)
149    REAL(rstd),INTENT(IN) :: ps(iim*jjm)
150    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm)
151    REAL(rstd),INTENT(OUT):: dtheta_rhodz(iim*jjm,llm)
152    REAL(rstd),INTENT(OUT):: dps(iim*jjm)
153   
154    INTEGER :: i,j,ij,l
155    REAL(rstd) :: ww,uu 
156    REAL(rstd) :: delta
157    REAL(rstd) :: etav,hv   
158    REAL(rstd),ALLOCATABLE,SAVE :: theta(:,:)  ! potential temperature
159    REAL(rstd),ALLOCATABLE,SAVE :: p(:,:)  ! pression
160    REAL(rstd),ALLOCATABLE,SAVE :: pk(:,:)   ! Exner function
161    REAL(rstd),ALLOCATABLE,SAVE :: pks(:)
162 ! Intermediate variable to compute exner function
163    REAL(rstd),ALLOCATABLE,SAVE :: alpha(:,:)
164    REAL(rstd),ALLOCATABLE,SAVE :: beta(:,:)
165 !   
166    REAL(rstd),ALLOCATABLE,SAVE :: phi(:,:)    ! geopotential
167    REAL(rstd),ALLOCATABLE,SAVE :: mass(:,:)   ! mass   
168    REAL(rstd),ALLOCATABLE,SAVE :: rhodz(:,:)   ! mass density   
169    REAL(rstd),ALLOCATABLE,SAVE :: Fe(:,:)   ! mass flux   
170    REAL(rstd),ALLOCATABLE,SAVE :: Ftheta(:,:) ! theta flux   
171    REAL(rstd),ALLOCATABLE,SAVE :: convm(:,:)    ! mass flux convergence
172    REAL(rstd),ALLOCATABLE,SAVE :: w(:,:)       ! vertical velocity     
173    REAL(rstd),ALLOCATABLE,SAVE :: qv(:,:)       ! potential velocity 
174    REAL(rstd),ALLOCATABLE,SAVE :: berni(:,:)   ! bernouilli term
175   
176    LOGICAL,SAVE :: first=.TRUE.
177!$OMP THREADPRIVATE(first)
178       
179!$OMP BARRIER     
180!$OMP MASTER 
181    IF (first) THEN
182      ALLOCATE(theta(iim*jjm,llm))  ! potential temperature
183      ALLOCATE(p(iim*jjm,llm+1))  ! pression
184      ALLOCATE(pk(iim*jjm,llm))   ! Exner function
185      ALLOCATE(pks(iim*jjm))
186      ALLOCATE(alpha(iim*jjm,llm))
187      ALLOCATE(beta(iim*jjm,llm))
188      ALLOCATE(phi(iim*jjm,llm))    ! geopotential
189      ALLOCATE(mass(iim*jjm,llm))   ! mass   
190      ALLOCATE(rhodz(iim*jjm,llm))   ! mass density   
191      ALLOCATE(Fe(3*iim*jjm,llm))   ! mass flux   
192      ALLOCATE(Ftheta(3*iim*jjm,llm)) ! theta flux   
193      ALLOCATE(convm(iim*jjm,llm))    ! mass flux convergence
194      ALLOCATE(w(iim*jjm,llm))       ! vertical velocity     
195      ALLOCATE(qv(2*iim*jjm,llm))       ! potential velocity 
196      ALLOCATE(berni(iim*jjm,llm))   ! bernouilli term
197      first=.FALSE.
198    ENDIF
199!$OMP END MASTER
200!$OMP BARRIER       
201   
202 !!! Compute pression
203   dtheta_rhodz(:,:)=0.
204   du(:,:)=0.
205 
206    DO    l    = 1, llm+1
207!$OMP DO
208      DO j=jj_begin-1,jj_end+1
209        DO i=ii_begin-1,ii_end+1
210          ij=(j-1)*iim+i
211          p(ij,l) = ap(l) + bp(l) * ps(ij)
212        ENDDO
213      ENDDO
214    ENDDO
215     
216 
217!!! Compute mass
218   DO l = 1, llm
219!$OMP DO
220     DO j=jj_begin-1,jj_end+1
221       DO i=ii_begin-1,ii_end+1
222         ij=(j-1)*iim+i
223         mass(ij,l) = ( p(ij,l) - p(ij,l+1) ) * Ai(ij)/g
224         rhodz(ij,l) = mass(ij,l) / Ai(ij)
225       ENDDO
226     ENDDO
227   ENDDO
228
229   
230!!!  Compute mass flux
231!! question à thomas : meilleure pondération de la masse sur les liens ?
232
233  DO l = 1, llm
234!$OMP DO
235    DO j=jj_begin-1,jj_end+1
236      DO i=ii_begin-1,ii_end+1
237        ij=(j-1)*iim+i
238        Fe(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)
239        Fe(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)
240        Fe(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown)
241      ENDDO
242    ENDDO
243  ENDDO
244
245
246!!! mass flux convergence computation 
247
248 ! horizontal convergence
249  DO l = 1, llm
250!$OMP DO
251    DO j=jj_begin,jj_end
252      DO i=ii_begin,ii_end
253        ij=(j-1)*iim+i
254!signe ?
255        convm(ij,l)= 1./Ai(ij)*(ne(ij,right)*Fe(ij+u_right,l)   +  &
256                                ne(ij,rup)*Fe(ij+u_rup,l)       +  & 
257                                ne(ij,lup)*Fe(ij+u_lup,l)       +  & 
258                                ne(ij,left)*Fe(ij+u_left,l)     +  & 
259                                ne(ij,ldown)*Fe(ij+u_ldown,l)   +  & 
260                                ne(ij,rdown)*Fe(ij+u_rdown,l))   
261       ENDDO
262     ENDDO
263   ENDDO
264 
265   
266 ! vertical integration from up to down
267  DO  l = llm-1, 1, -1
268!$OMP DO
269    DO j=jj_begin,jj_end
270      DO i=ii_begin,ii_end
271        ij=(j-1)*iim+i
272        convm(ij,l) = convm(ij,l) + convm(ij,l+1)
273      ENDDO
274    ENDDO
275  ENDDO
276 
277
278!!! Compute dps
279!$OMP DO
280  DO j=jj_begin,jj_end
281    DO i=ii_begin,ii_end
282      ij=(j-1)*iim+i
283      dps(ij)=-convm(ij,1) * g 
284    ENDDO
285  ENDDO
286                           
287  END SUBROUTINE compute_caldyn
288 
289 
290 
291END MODULE caldyn_adv_mod
Note: See TracBrowser for help on using the repository browser.