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

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

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