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

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