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

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

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