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

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

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