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

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

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