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

Last change on this file since 133 was 129, checked in by dubos, 11 years ago

Changed timeloop design for multistage schemes
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(write_out,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  LOGICAL,INTENT(IN)    :: write_out
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 (write_out) 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.