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

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

new hydrostatic balance - tested with JW06 - code cleanup to follow
File size: 34.1 KB
Line 
1MODULE caldyn_gcm_mod
2  USE icosa
3  USE transfert_mod
4  PRIVATE
5
6  INTEGER, PARAMETER :: energy=1, enstrophy=2
7  TYPE(t_field),POINTER :: f_out_u(:), f_p(:), f_rhodz(:), f_qu(:)
8  REAL(rstd),POINTER :: out_u(:,:), p(:,:), rhodz(:,:), qu(:,:)
9
10  TYPE(t_field),POINTER :: f_buf_i(:), f_buf_ulon(:), f_buf_ulat(:), f_buf_u3d(:)
11  TYPE(t_field),POINTER :: f_buf_v(:), f_buf_s(:), f_buf_p(:)
12
13! temporary shared variable for caldyn
14  TYPE(t_field),POINTER :: f_theta(:)
15  TYPE(t_field),POINTER :: f_pk(:)
16!  TYPE(t_field),POINTER :: f_pks(:)
17  TYPE(t_field),POINTER :: f_phi(:)
18  TYPE(t_field),POINTER :: f_geopot(:)
19  TYPE(t_field),POINTER :: f_divm(:)
20  TYPE(t_field),POINTER :: f_wwuu(:)
21   
22  INTEGER :: caldyn_hydrostat, caldyn_conserv
23 
24  TYPE(t_message) :: req_ps, req_theta_rhodz, req_u, req_qu
25
26  PUBLIC init_caldyn, caldyn, write_output_fields,req_ps
27
28CONTAINS
29 
30  SUBROUTINE init_caldyn
31    USE icosa
32    USE exner_mod
33    USE mpipara
34    IMPLICIT NONE
35    CHARACTER(len=255) :: def
36 
37    def='enstrophy'
38    CALL getin('caldyn_conserv',def)
39    SELECT CASE(TRIM(def))
40    CASE('energy')
41       caldyn_conserv=energy
42    CASE('enstrophy')
43       caldyn_conserv=enstrophy
44    CASE DEFAULT
45       IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_conserv : <', TRIM(def),'> options are <energy>, <enstrophy>'
46       STOP
47    END SELECT
48    IF (is_mpi_root) PRINT *, 'caldyn_conserv=',def
49
50    def='direct'
51    CALL getin('caldyn_exner',def)
52    SELECT CASE(TRIM(def))
53    CASE('lmdz')
54       caldyn_exner=lmdz
55    CASE('direct')
56       caldyn_exner=direct
57    CASE DEFAULT
58       IF (is_mpi_root) PRINT*,'Bad selector for variable caldyn_exner : <', TRIM(def),'> options are <lmdz>, <direct>'
59       STOP
60    END SELECT
61
62    def='direct'
63    CALL getin('caldyn_hydrostat',def)
64    SELECT CASE(TRIM(def))
65    CASE('lmdz')
66       caldyn_hydrostat=lmdz
67    CASE('direct')
68       caldyn_hydrostat=direct
69    CASE DEFAULT
70       IF (is_mpi_root) PRINT*,'Bad selector for variable caldyn_hydrostat : <', TRIM(def),'> options are <lmdz>, <direct>'
71       STOP
72    END SELECT
73   
74    CALL allocate_caldyn
75 
76  END SUBROUTINE init_caldyn
77 
78  SUBROUTINE allocate_caldyn
79  USE icosa
80  IMPLICIT NONE
81
82    CALL allocate_field(f_out_u,field_u,type_real,llm) 
83    CALL allocate_field(f_p,field_t,type_real,llm+1)
84    CALL allocate_field(f_rhodz,field_t,type_real,llm) 
85    CALL allocate_field(f_qu,field_u,type_real,llm) 
86 
87    CALL allocate_field(f_buf_i,field_t,type_real,llm)
88    CALL allocate_field(f_buf_p,field_t,type_real,llm+1) 
89    CALL allocate_field(f_buf_u3d,field_t,type_real,3,llm)  ! 3D vel at cell centers
90    CALL allocate_field(f_buf_ulon,field_t,type_real,llm)
91    CALL allocate_field(f_buf_ulat,field_t,type_real,llm)
92    CALL allocate_field(f_buf_v,field_z,type_real,llm)
93    CALL allocate_field(f_buf_s,field_t,type_real)
94
95    CALL allocate_field(f_theta,field_t,type_real,llm)  ! potential temperature
96    CALL allocate_field(f_pk,field_t,type_real,llm)
97!    CALL allocate_field(f_pks,field_t,type_real) ! Exner function
98!    CALL allocate_field(f_phi,field_t,type_real,llm)     ! geopotential
99    CALL allocate_field(f_geopot,field_t,type_real,llm+1)  ! geopotential (new)
100    CALL allocate_field(f_divm,field_t,type_real,llm)  ! mass flux divergence
101    CALL allocate_field(f_wwuu,field_u,type_real,llm+1) 
102
103  END SUBROUTINE allocate_caldyn
104   
105  SUBROUTINE caldyn(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, &
106       f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du)
107    USE icosa
108    USE vorticity_mod
109    USE kinetic_mod
110    USE theta2theta_rhodz_mod
111    USE mpipara
112    USE trace
113    USE omp_para
114    IMPLICIT NONE
115    LOGICAL,INTENT(IN)    :: write_out
116    TYPE(t_field),POINTER :: f_phis(:)
117    TYPE(t_field),POINTER :: f_ps(:)
118    TYPE(t_field),POINTER :: f_theta_rhodz(:)
119    TYPE(t_field),POINTER :: f_u(:)
120    TYPE(t_field),POINTER :: f_q(:)
121    TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:)
122    TYPE(t_field),POINTER :: f_dps(:)
123    TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
124    TYPE(t_field),POINTER :: f_du(:)
125   
126    REAL(rstd),POINTER :: phis(:), ps(:), dps(:)
127    REAL(rstd),POINTER :: theta_rhodz(:,:), dtheta_rhodz(:,:)
128    REAL(rstd),POINTER :: u(:,:), du(:,:), hflux(:,:), wflux(:,:)
129    REAL(rstd),POINTER :: p(:,:), rhodz(:,:), qu(:,:)
130
131! temporary shared variable
132    REAL(rstd),POINTER  :: theta(:,:) 
133    REAL(rstd),POINTER  :: pk(:,:) !, pks(:)
134    REAL(rstd),POINTER  :: phi(:,:)   
135    REAL(rstd),POINTER  :: geopot(:,:)
136    REAL(rstd),POINTER  :: divm(:,:) 
137    REAL(rstd),POINTER  :: wwuu(:,:)
138   
139   
140    INTEGER :: ind,ij
141    LOGICAL,SAVE :: first=.TRUE.
142!$OMP THREADPRIVATE(first)
143
144   
145    IF (first) THEN
146      first=.FALSE.
147      CALL init_message(f_ps,req_i1,req_ps)
148      CALL init_message(f_theta_rhodz,req_i1,req_theta_rhodz)
149      CALL init_message(f_u,req_e1_vect,req_u)
150      CALL init_message(f_qu,req_e1_scal,req_qu)
151      CALL send_message(f_ps,req_ps) 
152    ENDIF
153   
154    CALL trace_start("caldyn")
155
156    CALL send_message(f_u,req_u)
157    CALL send_message(f_theta_rhodz,req_theta_rhodz) 
158
159    SELECT CASE(caldyn_conserv)
160    CASE(energy) ! energy-conserving
161       DO ind=1,ndomain
162          CALL swap_dimensions(ind)
163          CALL swap_geometry(ind)
164          ps=f_ps(ind)
165          rhodz=f_rhodz(ind)
166          p=f_p(ind)
167          qu=f_qu(ind)
168          u=f_u(ind)
169
170          CALL compute_pvort(ps, u, p,rhodz,qu)
171
172       ENDDO
173
174       CALL send_message(f_qu,req_qu)
175
176       DO ind=1,ndomain
177          CALL swap_dimensions(ind)
178          CALL swap_geometry(ind)
179          phis=f_phis(ind)
180          hflux=f_hflux(ind)
181          wflux=f_wflux(ind)
182          ps=f_ps(ind)
183          dps=f_dps(ind)
184          theta_rhodz=f_theta_rhodz(ind)
185          dtheta_rhodz=f_dtheta_rhodz(ind)
186          rhodz=f_rhodz(ind)
187          p=f_p(ind)
188          qu=f_qu(ind)
189          u=f_u(ind)
190          du=f_du(ind)
191          out_u=f_out_u(ind)
192          theta = f_theta(ind)
193          pk = f_pk(ind)
194!          pks = f_pks(ind)
195!          phi = f_phi(ind)
196          geopot = f_geopot(ind) 
197          divm = f_divm(ind)
198          wwuu=f_wwuu(ind)
199
200          CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, hflux, wflux, dps, dtheta_rhodz, du, &
201                              theta,pk, geopot, divm, wwuu)
202!                              theta,pk, pks, phi, geopot, divm, wwuu)
203       ENDDO       
204       
205    CASE(enstrophy) ! enstrophy-conserving
206       DO ind=1,ndomain
207          CALL swap_dimensions(ind)
208          CALL swap_geometry(ind)
209          phis=f_phis(ind)
210          phi=f_phi(ind)
211          geopot = f_geopot(ind) 
212          ps=f_ps(ind)
213          dps=f_dps(ind)
214          hflux=f_hflux(ind)
215          wflux=f_wflux(ind)
216          theta_rhodz=f_theta_rhodz(ind)
217          dtheta_rhodz=f_dtheta_rhodz(ind)
218          rhodz=f_rhodz(ind)
219          p=f_p(ind)
220          qu=f_qu(ind)
221          u=f_u(ind)
222          du=f_du(ind)
223          out_u=f_out_u(ind) 
224          wwuu=f_wwuu(ind)
225          CALL compute_pvort(ps, u, p,rhodz,qu)
226          CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, hflux, wflux, dps, dtheta_rhodz, du,   &
227                              theta,pk, phi, divm, wwuu)
228!                              theta,pk, pks, phi, geopot, divm, wwuu)
229       ENDDO
230       
231    CASE DEFAULT
232       STOP
233    END SELECT
234
235!$OMP BARRIER
236    IF (write_out) THEN
237
238       IF (is_mpi_root) PRINT *,'CALL write_output_fields'
239
240! ---> for openMP test to fix later
241!       CALL write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
242!            f_buf_i, f_buf_v, f_buf_u3d, f_buf_ulon, f_buf_ulat, f_buf_s, f_buf_p)
243       CALL writefield("ps",f_ps)
244       CALL writefield("dps",f_dps)
245       CALL writefield("vort",f_qu)
246       CALL writefield("theta",f_theta)
247
248    END IF
249   
250    !    CALL check_mass_conservation(f_ps,f_dps)
251    CALL trace_end("caldyn")
252!$OMP BARRIER
253   
254END SUBROUTINE caldyn
255   
256SUBROUTINE compute_pvort(ps, u, p,rhodz,qu)
257  USE icosa
258  USE disvert_mod
259  USE exner_mod
260  USE trace
261  USE omp_para
262  IMPLICIT NONE
263  REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
264  REAL(rstd),INTENT(IN)  :: ps(iim*jjm)
265  REAL(rstd),INTENT(OUT) :: p(iim*jjm,llm+1)
266  REAL(rstd),INTENT(OUT) :: rhodz(iim*jjm,llm)
267  REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm)
268 
269  INTEGER :: i,j,ij,l
270  REAL(rstd) :: etav,hv
271  REAL(rstd) :: qv(2*iim*jjm,llm)     ! potential velocity 
272 
273
274  CALL trace_start("compute_pvort") 
275
276  CALL wait_message(req_ps)   
277!!! Compute pressure
278   DO    l    = ll_begin, ll_endp1
279    CALL test_message(req_u) 
280
281     DO j=jj_begin-1,jj_end+1
282        DO i=ii_begin-1,ii_end+1
283           ij=(j-1)*iim+i
284           p(ij,l) = ap(l) + bp(l) * ps(ij)
285        ENDDO
286     ENDDO
287  ENDDO
288
289!$OMP BARRIER 
290
291!!! Compute mass
292   DO l = ll_begin,ll_end
293     CALL test_message(req_u) 
294     DO j=jj_begin-1,jj_end+1
295        DO i=ii_begin-1,ii_end+1
296           ij=(j-1)*iim+i
297           rhodz(ij,l) = ( p(ij,l) - p(ij,l+1) )/g
298        ENDDO
299     ENDDO
300  ENDDO
301
302  CALL wait_message(req_u)   
303 
304!!! Compute shallow-water potential vorticity
305  DO l = ll_begin,ll_end
306    CALL test_message(req_theta_rhodz) 
307
308    DO j=jj_begin-1,jj_end+1
309       DO i=ii_begin-1,ii_end+1
310          ij=(j-1)*iim+i
311         
312          etav= 1./Av(ij+z_up)*(  ne_rup        * u(ij+u_rup,l)        * de(ij+u_rup)         &
313                                + ne_left * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  &
314                                - ne_lup        * u(ij+u_lup,l)        * de(ij+u_lup) )                               
315
316          hv =  Riv2(ij,vup)          * rhodz(ij,l)            &
317              + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
318              + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
319     
320          qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv
321         
322          etav = 1./Av(ij+z_down)*(  ne_ldown         * u(ij+u_ldown,l)          * de(ij+u_ldown)          &
323                                   + ne_right * u(ij+t_ldown+u_right,l)  * de(ij+t_ldown+u_right)  &
324                                   - ne_rdown         * u(ij+u_rdown,l)          * de(ij+u_rdown) )
325       
326          hv =  Riv2(ij,vdown)        * rhodz(ij,l)          &
327              + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
328              + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
329                       
330          qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv
331
332        ENDDO
333      ENDDO
334
335      DO j=jj_begin,jj_end
336         DO i=ii_begin,ii_end
337            ij=(j-1)*iim+i
338            qu(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l)) 
339            qu(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l)) 
340            qu(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l)) 
341         END DO
342      END DO
343     
344   ENDDO
345
346   CALL trace_end("compute_pvort")
347                                                       
348  END SUBROUTINE compute_pvort
349
350 
351  SUBROUTINE compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, hflux, wflux, dps, dtheta_rhodz, du,   &
352                            theta,pk, geopot, divm, wwuu) 
353!                            theta,pk, pks, phi, geopot, divm, wwuu)
354  USE icosa
355  USE disvert_mod
356  USE exner_mod
357  USE trace
358  USE omp_para
359  IMPLICIT NONE
360    REAL(rstd),INTENT(IN)  :: phis(iim*jjm)
361    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
362    REAL(rstd),INTENT(IN)  :: theta_rhodz(iim*jjm,llm)
363    REAL(rstd),INTENT(IN)  :: ps(iim*jjm)
364    REAL(rstd),INTENT(IN)  :: p(iim*jjm,llm+1)
365    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
366    REAL(rstd),INTENT(IN)  :: qu(iim*3*jjm,llm)
367
368    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm), hflux(iim*3*jjm,llm) ! hflux in kg/s
369    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm)
370    REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
371    REAL(rstd),INTENT(OUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s)
372
373! temporary variable
374    REAL(rstd),INTENT(INOUT) :: theta(iim*jjm,llm)  ! potential temperature
375    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) !, pks(iim*jjm) ! Exner function
376!    REAL(rstd),INTENT(INOUT) :: phi(iim*jjm,llm)    ! geopotential
377    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)    ! geopotential
378    REAL(rstd),INTENT(INOUT) :: divm(iim*jjm,llm)  ! mass flux divergence
379    REAL(rstd),INTENT(INOUT) :: wwuu(iim*3*jjm,llm+1)
380   
381    REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux
382    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernouilli function
383   
384    INTEGER :: i,j,ij,l
385    REAL(rstd) :: ww,uu, p_ik, exner_ik
386
387    CALL trace_start("compute_caldyn")
388
389    CALL wait_message(req_theta_rhodz) 
390
391!!! Compute theta
392   DO l = ll_begin, ll_end
393      IF (caldyn_conserv==energy) CALL test_message(req_qu) 
394      DO j=jj_begin-1,jj_end+1
395         DO i=ii_begin-1,ii_end+1
396            ij=(j-1)*iim+i
397            theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
398         ENDDO
399      ENDDO
400   ENDDO
401
402  DO l = ll_begin, ll_end
403!!!  Compute mass and theta fluxes
404    IF (caldyn_conserv==energy) CALL test_message(req_qu) 
405    DO j=jj_begin-1,jj_end+1
406      DO i=ii_begin-1,ii_end+1
407        ij=(j-1)*iim+i
408        hflux(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)
409        hflux(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)
410        hflux(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown)
411
412        Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*hflux(ij+u_right,l)
413        Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*hflux(ij+u_lup,l)
414        Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*hflux(ij+u_ldown,l)
415      ENDDO
416    ENDDO
417
418!!! compute horizontal divergence of fluxes
419    DO j=jj_begin,jj_end
420      DO i=ii_begin,ii_end
421        ij=(j-1)*iim+i 
422        ! divm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
423        divm(ij,l)= 1./Ai(ij)*(ne_right*hflux(ij+u_right,l)   +  &
424                                ne_rup*hflux(ij+u_rup,l)       +  & 
425                                ne_lup*hflux(ij+u_lup,l)       +  & 
426                                ne_left*hflux(ij+u_left,l)     +  & 
427                                ne_ldown*hflux(ij+u_ldown,l)   +  & 
428                                ne_rdown*hflux(ij+u_rdown,l))   
429
430        ! signe ? attention d (rho theta dz)
431        ! dtheta_rhodz =  -div(flux.theta)
432        dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne_right*Ftheta(ij+u_right,l)    +  &
433                                 ne_rup*Ftheta(ij+u_rup,l)        +  & 
434                                 ne_lup*Ftheta(ij+u_lup,l)        +  & 
435                                 ne_left*Ftheta(ij+u_left,l)      +  & 
436                                 ne_ldown*Ftheta(ij+u_ldown,l)    +  & 
437                                 ne_rdown*Ftheta(ij+u_rdown,l))
438      ENDDO
439    ENDDO
440  ENDDO
441
442!$OMP BARRIER   
443!!! cumulate mass flux divergence from top to bottom
444  DO  l = llm-1, 1, -1
445    IF (caldyn_conserv==energy) CALL test_message(req_qu) 
446!$OMP DO SCHEDULE(STATIC)
447    DO j=jj_begin,jj_end
448      DO i=ii_begin,ii_end
449        ij=(j-1)*iim+i
450        divm(ij,l) = divm(ij,l) + divm(ij,l+1)
451      ENDDO
452    ENDDO
453  ENDDO
454
455! IMPLICIT FLUSH on divm
456!!!!!!!!!!!!!!!!!!!!!!!!!
457 
458!!! Compute vertical mass flux
459!  DO l = 2,llm
460  DO l=ll_beginp1,ll_end
461    IF (caldyn_conserv==energy) CALL test_message(req_qu) 
462    DO j=jj_begin,jj_end
463      DO i=ii_begin,ii_end
464        ij=(j-1)*iim+i
465        ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
466        ! => w>0 for upward transport
467        wflux( ij, l ) = divm( ij, l ) - bp(l) * divm( ij, 1 )
468      ENDDO
469    ENDDO
470  ENDDO
471
472! compute dps, set vertical mass flux at the surface to 0
473  IF (omp_first) THEN
474    DO j=jj_begin,jj_end
475      DO i=ii_begin,ii_end
476        ij=(j-1)*iim+i
477        wflux(ij,1)  = 0.
478        ! dps/dt = -int(div flux)dz
479        dps(ij)=-divm(ij,1) * g 
480      ENDDO
481    ENDDO
482  ENDIF
483
484
485  IF (omp_last) THEN
486    DO j=jj_begin,jj_end
487      DO i=ii_begin,ii_end
488        ij=(j-1)*iim+i
489        wflux(ij,llm+1)  = 0.
490      ENDDO
491    ENDDO
492  ENDIF
493!!! Compute potential vorticity (Coriolis) contribution to du
494
495    SELECT CASE(caldyn_conserv)
496    CASE(energy) ! energy-conserving TRiSK
497
498       CALL wait_message(req_qu)
499
500        DO l=ll_begin,ll_end
501          DO j=jj_begin,jj_end
502             DO i=ii_begin,ii_end
503                ij=(j-1)*iim+i
504
505                uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             &
506                     wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             &
507                     wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+                           &
508                     wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+                         &
509                     wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+                         & 
510                     wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_ldown,l))+         &
511                     wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rdown,l))+         &
512                     wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))+         &
513                     wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rup,l))+             &
514                     wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_lup,l))
515                du(ij+u_right,l) = .5*uu/de(ij+u_right)
516               
517                uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        &
518                     wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) +                      &
519                     wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +                      &
520                     wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) +                      &
521                     wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) +                          & 
522                     wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_right,l)) +          &
523                     wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_rup,l)) +              &
524                     wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) +              &
525                     wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_left,l)) +            &
526                     wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_ldown,l))
527                du(ij+u_lup,l) = .5*uu/de(ij+u_lup)
528
529               
530                uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      &
531                     wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) +                      &
532                     wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +                          &
533                     wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) +                          &
534                     wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) +                        & 
535                     wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_lup,l)) +          &
536                     wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_left,l)) +        &
537                     wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) +      &
538                     wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l)) +      &
539                     wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_right,l))
540                du(ij+u_ldown,l) = .5*uu/de(ij+u_ldown)
541               
542             ENDDO
543          ENDDO
544       ENDDO
545       
546    CASE(enstrophy) ! enstrophy-conserving TRiSK
547 
548        DO l=ll_begin,ll_end
549          DO j=jj_begin,jj_end
550             DO i=ii_begin,ii_end
551                ij=(j-1)*iim+i
552
553                uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+                           &
554                     wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+                           &
555                     wee(ij+u_right,3,1)*hflux(ij+u_left,l)+                          &
556                     wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+                         &
557                     wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+                         & 
558                     wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+                 &
559                     wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+                 &
560                     wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+                 &
561                     wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+                   &
562                     wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)
563                du(ij+u_right,l) = qu(ij+u_right,l)*uu/de(ij+u_right)
564               
565               
566                uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+                        &
567                     wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+                       &
568                     wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+                       &
569                     wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+                       &
570                     wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+                         & 
571                     wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+                 &
572                     wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+                   &
573                     wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+                   &
574                     wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+                  &
575                     wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)
576                du(ij+u_lup,l) = qu(ij+u_lup,l)*uu/de(ij+u_lup)
577
578                uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+                      &
579                     wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+                      &
580                     wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+                        &
581                     wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+                        &
582                     wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+                       & 
583                     wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+                &
584                     wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+               &
585                     wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+              &
586                     wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+              &
587                     wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)
588                du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu/de(ij+u_ldown)
589     
590             ENDDO
591          ENDDO
592       ENDDO
593       
594    CASE DEFAULT
595       STOP
596    END SELECT
597 
598!!!! Compute Exner function
599!!    CALL compute_exner(ps,p,pks,pk,1)
600!! replaced in source
601!    IF (omp_first) THEN   
602!       DO j=jj_begin-1,jj_end+1
603!          DO i=ii_begin-1,ii_end+1
604!             ij=(j-1)*iim+i
605!             pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
606!          ENDDO
607!       ENDDO
608!     ENDIF
609!     
610!       ! 3D : pk
611!     DO l = ll_begin, ll_end
612!        DO j=jj_begin-1,jj_end+1
613!           DO i=ii_begin-1,ii_end+1
614!              ij=(j-1)*iim+i
615!              pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa
616!           ENDDO
617!        ENDDO
618!     ENDDO
619!
620!!---> flush pk,pks, theta
621!!$OMP BARRIER
622!
623!!! Compute geopotential (old)
624!
625!  ! for first layer
626!  IF (omp_first) THEN
627!    DO j=jj_begin-1,jj_end+1
628!      DO i=ii_begin-1,ii_end+1
629!        ij=(j-1)*iim+i
630!        phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
631!      ENDDO
632!    ENDDO
633!  ENDIF
634!!!-> implicit flush on phi(:,1)
635!         
636!  ! for other layers
637!   DO l = ll_beginp1, ll_end
638!     DO j=jj_begin-1,jj_end+1
639!       DO i=ii_begin-1,ii_end+1
640!         ij=(j-1)*iim+i
641!         phi(ij,l) =  0.5 * ( theta(ij,l)  + theta(ij,l-1) )  &
642!                          * (  pk(ij,l-1) -  pk(ij,l)    )
643!       ENDDO
644!     ENDDO
645!   ENDDO
646!
647!!$OMP BARRIER
648!   DO l = 2, llm
649!!$OMP DO
650!     DO j=jj_begin-1,jj_end+1
651!       DO i=ii_begin-1,ii_end+1
652!         ij=(j-1)*iim+i
653!         phi(ij,l) = phi(ij,l)+ phi(ij,l-1)
654!       ENDDO
655!     ENDDO
656!   ENDDO
657!! --> IMPLICIT FLUSH on phi
658     
659!!! Compute exner function and geopotential
660
661! geopot=phis for first layer
662!$OMP DO SCHEDULE(STATIC)
663    DO j=jj_begin-1,jj_end+1
664       DO i=ii_begin-1,ii_end+1
665          ij=(j-1)*iim+i
666          geopot(ij,1) = phis(ij)
667       ENDDO
668    ENDDO
669! for other layers
670    DO l = 1,llm
671!$OMP DO SCHEDULE(STATIC)
672       DO j=jj_begin-1,jj_end+1
673          DO i=ii_begin-1,ii_end+1
674             ij=(j-1)*iim+i
675             p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later
676             !         p_ik = ptop + g*(mass_ak(l)+ mass_bk(l)*ps(i,j))
677             exner_ik = cpp * (p_ik/preff) ** kappa
678             pk(ij,l) = exner_ik
679             ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
680             geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 
681          ENDDO
682       ENDDO
683    ENDDO
684
685!!! Compute bernouilli term = Kinetic Energy + geopotential
686   DO l=ll_begin,ll_end
687    DO j=jj_begin,jj_end
688      DO i=ii_begin,ii_end
689        ij=(j-1)*iim+i
690
691        berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
692                     + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
693                                     le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
694                                     le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
695                                     le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       &
696                                     le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
697                                     le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 
698       
699      ENDDO
700    ENDDO
701  ENDDO
702
703 
704!!! gradients of Bernoulli and Exner functions
705
706   DO l=ll_begin,ll_end
707    DO j=jj_begin,jj_end
708      DO i=ii_begin,ii_end
709        ij=(j-1)*iim+i
710       
711        du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right) * (             &
712                          0.5*(theta(ij,l)+theta(ij+t_right,l))                &
713                             *( ne_right*pk(ij,l)+ne_left*pk(ij+t_right,l))    &
714                        + ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l) )
715       
716       
717        du(ij+u_lup,l) = du(ij+u_lup,l) +  1/de(ij+u_lup) * (                  &
718                          0.5*(theta(ij,l)+theta(ij+t_lup,l))                  &
719                             *( ne_lup*pk(ij,l)+ne_rdown*pk(ij+t_lup,l))       &
720                        + ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) )
721               
722        du(ij+u_ldown,l) = du(ij+u_ldown,l) + 1/de(ij+u_ldown) * (             &
723                          0.5*(theta(ij,l)+theta(ij+t_ldown,l))                &
724                             *( ne_ldown*pk(ij,l)+ne_rup*pk(ij+t_ldown,l))     &
725                        + ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) )
726
727      ENDDO
728    ENDDO
729  ENDDO
730
731   
732  DO l=ll_begin,ll_endm1
733    DO j=jj_begin,jj_end
734      DO i=ii_begin,ii_end
735        ij=(j-1)*iim+i
736        dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -   0.5 * (  wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1)))
737      ENDDO
738    ENDDO
739  ENDDO
740 
741  DO l=ll_beginp1,ll_end
742    DO j=jj_begin,jj_end
743      DO i=ii_begin,ii_end
744        ij=(j-1)*iim+i
745        dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  +   0.5 * ( wflux(ij,l  ) * (theta(ij,l-1) + theta(ij,l) ) )
746      ENDDO
747    ENDDO
748  ENDDO
749 
750  IF (omp_first) THEN
751    DO j=jj_begin,jj_end
752      DO i=ii_begin,ii_end
753        ij=(j-1)*iim+i
754        wwuu(ij+u_right,1)=0   
755        wwuu(ij+u_lup,1)=0   
756        wwuu(ij+u_ldown,1)=0
757      ENDDO
758    ENDDO
759  ENDIF   
760
761  IF (omp_last) THEN
762    DO j=jj_begin,jj_end
763      DO i=ii_begin,ii_end
764        ij=(j-1)*iim+i
765        wwuu(ij+u_right,llm+1)=0   
766        wwuu(ij+u_lup,llm+1)=0   
767        wwuu(ij+u_ldown,llm+1)=0
768      ENDDO
769    ENDDO
770  ENDIF   
771 
772  DO l=ll_beginp1,ll_end
773    DO j=jj_begin,jj_end
774      DO i=ii_begin,ii_end
775        ij=(j-1)*iim+i
776        wwuu(ij+u_right,l) = 0.5*( wflux(ij,l) + wflux(ij+t_right,l)) * (u(ij+u_right,l) - u(ij+u_right,l-1))
777        wwuu(ij+u_lup,l) = 0.5* ( wflux(ij,l) + wflux(ij+t_lup,l)) * (u(ij+u_lup,l) - u(ij+u_lup,l-1))
778        wwuu(ij+u_ldown,l) = 0.5*( wflux(ij,l) + wflux(ij+t_ldown,l)) * (u(ij+u_ldown,l) - u(ij+u_ldown,l-1))
779       ENDDO
780     ENDDO
781   ENDDO
782
783 !--> flush wwuu 
784 !$OMP BARRIER
785
786  DO l=ll_begin,ll_end
787    DO j=jj_begin,jj_end
788      DO i=ii_begin,ii_end
789        ij=(j-1)*iim+i
790        du(ij+u_right, l )   = du(ij+u_right,l)  - (wwuu(ij+u_right,l+1)+ wwuu(ij+u_right,l)) / (rhodz(ij,l)+rhodz(ij+t_right,l))
791        du(ij+u_lup, l )     = du(ij+u_lup,l)    - (wwuu(ij+u_lup,l+1)  + wwuu(ij+u_lup,l))   / (rhodz(ij,l)+rhodz(ij+t_lup,l))
792        du(ij+u_ldown, l )   = du(ij+u_ldown,l)  - (wwuu(ij+u_ldown,l+1)+ wwuu(ij+u_ldown,l)) / (rhodz(ij,l)+rhodz(ij+t_ldown,l))
793      ENDDO
794    ENDDO
795  ENDDO     
796 
797  CALL trace_end("compute_caldyn")
798
799  END SUBROUTINE compute_caldyn
800
801!-------------------------------- Diagnostics ----------------------------
802
803  SUBROUTINE check_mass_conservation(f_ps,f_dps)
804  USE icosa
805  USE mpipara
806  IMPLICIT NONE
807    TYPE(t_field),POINTER :: f_ps(:)
808    TYPE(t_field),POINTER :: f_dps(:)
809    REAL(rstd),POINTER :: ps(:)
810    REAL(rstd),POINTER :: dps(:)
811    REAL(rstd) :: mass_tot,dmass_tot
812    INTEGER :: ind,i,j,ij
813   
814    mass_tot=0
815    dmass_tot=0
816   
817    CALL transfert_request(f_dps,req_i1)
818    CALL transfert_request(f_ps,req_i1)
819
820    DO ind=1,ndomain
821      CALL swap_dimensions(ind)
822      CALL swap_geometry(ind)
823
824      ps=f_ps(ind)
825      dps=f_dps(ind)
826
827      DO j=jj_begin,jj_end
828        DO i=ii_begin,ii_end
829          ij=(j-1)*iim+i
830          IF (domain(ind)%own(i,j)) THEN
831            mass_tot=mass_tot+ps(ij)*Ai(ij)/g
832            dmass_tot=dmass_tot+dps(ij)*Ai(ij)/g
833          ENDIF
834        ENDDO
835      ENDDO
836   
837    ENDDO
838    IF (is_mpi_root) PRINT*, "mass_tot ", mass_tot,"      dmass_tot ",dmass_tot       
839
840  END SUBROUTINE check_mass_conservation 
841 
842  SUBROUTINE write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
843       f_buf_i, f_buf_v, f_buf_i3, f_buf1_i, f_buf2_i, f_buf_s, f_buf_p)
844    USE icosa
845    USE vorticity_mod
846    USE theta2theta_rhodz_mod
847    USE pression_mod
848    USE omega_mod
849    USE write_field
850    USE vertical_interp_mod
851    USE wind_mod
852    TYPE(t_field),POINTER :: f_ps(:), f_phis(:), f_u(:), f_theta_rhodz(:), f_q(:), f_dps(:), &
853         f_buf_i(:), f_buf_v(:), f_buf_i3(:), f_buf1_i(:), f_buf2_i(:), f_buf_s(:), f_buf_p(:)
854   
855    REAL(rstd) :: out_pression_lev
856    CHARACTER(LEN=255) :: str_pression
857    CHARACTER(LEN=255) :: physics_type
858   
859    out_pression_level=0
860    CALL getin("out_pression_level",out_pression_level) 
861    WRITE(str_pression,*) INT(out_pression_level/100)
862    str_pression=ADJUSTL(str_pression)
863   
864    CALL writefield("ps",f_ps)
865    CALL writefield("dps",f_dps)
866    CALL writefield("phis",f_phis)
867    CALL vorticity(f_u,f_buf_v)
868    CALL writefield("vort",f_buf_v)
869
870    CALL w_omega(f_ps, f_u, f_buf_i)
871    CALL writefield('omega', f_buf_i)
872    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
873      CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
874      CALL writefield("omega"//TRIM(str_pression),f_buf_s)
875    ENDIF
876   
877    ! Temperature
878    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ;
879   
880    CALL getin('physics',physics_type)
881    IF (TRIM(physics_type)=='dcmip') THEN
882      CALL Tv2T(f_buf_i,f_q,f_buf1_i) 
883      CALL writefield("T",f_buf1_i)
884      IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
885        CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
886        CALL writefield("T"//TRIM(str_pression),f_buf_s)
887      ENDIF
888    ELSE
889      CALL writefield("T",f_buf_i)
890      IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
891        CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
892        CALL writefield("T"//TRIM(str_pression),f_buf_s)
893      ENDIF
894    ENDIF
895   
896    ! velocity components
897    CALL un2ulonlat(f_u, f_buf1_i, f_buf2_i)
898    CALL writefield("ulon",f_buf1_i)
899    CALL writefield("ulat",f_buf2_i)
900
901    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
902      CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
903      CALL writefield("ulon"//TRIM(str_pression),f_buf_s)
904      CALL vertical_interp(f_ps,f_buf2_i,f_buf_s,out_pression_level)
905      CALL writefield("ulat"//TRIM(str_pression),f_buf_s)
906    ENDIF
907   
908    ! geopotential
909    CALL thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_buf_s,f_buf_p,f_buf1_i,f_buf2_i,f_buf_i)
910    CALL writefield("p",f_buf_p)
911    CALL writefield("phi",f_buf_i)
912    CALL writefield("theta",f_buf1_i) ! potential temperature
913    CALL writefield("pk",f_buf2_i) ! Exner pressure
914 
915 
916  END SUBROUTINE write_output_fields
917 
918  SUBROUTINE thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_pks,f_p,f_theta,f_pk,f_phi) 
919    USE field_mod
920    USE pression_mod
921    USE exner_mod
922    USE geopotential_mod
923    USE theta2theta_rhodz_mod
924    TYPE(t_field), POINTER :: f_ps(:), f_phis(:), f_theta_rhodz(:), &  ! IN
925         f_pks(:), f_p(:), f_theta(:), f_pk(:), f_phi(:)               ! OUT
926    REAL(rstd),POINTER :: pk(:,:), p(:,:), theta(:,:), theta_rhodz(:,:), &
927         phi(:,:), phis(:), ps(:), pks(:)
928    INTEGER :: ind
929
930    DO ind=1,ndomain
931       CALL swap_dimensions(ind)
932       CALL swap_geometry(ind)
933       ps = f_ps(ind)
934       p  = f_p(ind)
935       CALL compute_pression(ps,p,0)
936       pk = f_pk(ind)
937       pks = f_pks(ind)
938       CALL compute_exner(ps,p,pks,pk,0)
939       theta_rhodz = f_theta_rhodz(ind)
940       theta = f_theta(ind)
941       CALL compute_theta_rhodz2theta(ps, theta_rhodz,theta,0)
942       phis = f_phis(ind)
943       phi = f_phi(ind)
944       CALL compute_geopotential(phis,pks,pk,theta,phi,0)
945    END DO
946
947  END SUBROUTINE thetarhodz2geopot
948 
949  SUBROUTINE Tv2T(f_Tv, f_q, f_T)
950  USE icosa
951  IMPLICIT NONE
952    TYPE(t_field), POINTER :: f_TV(:)
953    TYPE(t_field), POINTER :: f_q(:)
954    TYPE(t_field), POINTER :: f_T(:)
955   
956    REAL(rstd),POINTER :: Tv(:,:), q(:,:,:), T(:,:)
957    INTEGER :: ind
958   
959    DO ind=1,ndomain
960       CALL swap_dimensions(ind)
961       CALL swap_geometry(ind)
962       Tv=f_Tv(ind)
963       q=f_q(ind)
964       T=f_T(ind)
965       T=Tv/(1+0.608*q(:,:,1))
966    END DO
967   
968  END SUBROUTINE Tv2T
969 
970END MODULE caldyn_gcm_mod
Note: See TracBrowser for help on using the repository browser.