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

Last change on this file since 347 was 347, checked in by dubos, 9 years ago

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