source: trunk/ICOSA_LMDZ/src/phyvenus/interface_icosa_lmdz.f90 @ 3557

Last change on this file since 3557 was 3271, checked in by emillour, 9 months ago

Dynamico interfaces:
Minor fix to follow-up on recent dynamico updates.
"wind_mod" module no longer exist (but we didn't use it anyway).
EM

File size: 28.9 KB
Line 
1MODULE interface_icosa_lmdz_mod
2
3  USE field_mod, ONLY: t_field
4  USE transfert_mod, ONLY: t_message
5 
6 
7  TYPE(t_message),SAVE :: req_u
8  TYPE(t_message),SAVE :: req_dps0, req_dulon0, req_dulat0, req_dTemp0, req_dq0
9
10  TYPE(t_field),POINTER,SAVE :: f_p(:)
11  TYPE(t_field),POINTER,SAVE :: f_pks(:) 
12  TYPE(t_field),POINTER,SAVE :: f_pk(:) 
13  TYPE(t_field),POINTER,SAVE :: f_p_layer(:)   
14  TYPE(t_field),POINTER,SAVE :: f_theta(:)   
15  TYPE(t_field),POINTER,SAVE :: f_phi(:)   
16  TYPE(t_field),POINTER,SAVE :: f_Temp(:)   
17  TYPE(t_field),POINTER,SAVE :: f_ulon(:)   
18  TYPE(t_field),POINTER,SAVE :: f_ulat(:)   
19  TYPE(t_field),POINTER,SAVE :: f_dulon(:)
20  TYPE(t_field),POINTER,SAVE :: f_dulat(:)
21  TYPE(t_field),POINTER,SAVE :: f_dTemp(:)
22  TYPE(t_field),POINTER,SAVE :: f_dq(:)
23  TYPE(t_field),POINTER,SAVE :: f_dps(:)
24  TYPE(t_field),POINTER,SAVE :: f_duc(:)
25  TYPE(t_field),POINTER,SAVE :: f_bounds_lon(:)
26  TYPE(t_field),POINTER,SAVE :: f_bounds_lat(:)
27
28  INTEGER,SAVE :: start_clock
29  INTEGER,SAVE :: stop_clock
30  INTEGER,SAVE :: count_clock=0
31 
32  REAL,SAVE :: day_length ! length of a day (s)
33!  REAL,SAVE :: year_length ! length of a year (s)
34  INTEGER,SAVE :: start_day ! reference sol value at beginning of the run wrt Ls=0
35 
36  INTEGER,SAVE :: nbp_phys
37  INTEGER,SAVE :: nbp_phys_glo
38 
39  CHARACTER(len=30),SAVE,ALLOCATABLE :: tname(:) ! tracer names
40  CHARACTER(len=33),SAVE,ALLOCATABLE :: ttext(:) ! tracer long name for diagnostics
41  REAL,SAVE :: pday ! number of ellapsed sols since Ls=0
42  REAL,SAVE :: ptime ! "universal time" as fraction of sol (e.g. 0.5 for noon)
43
44
45CONTAINS
46
47  SUBROUTINE initialize_physics
48  USE distrib_icosa_lmdz_mod, ONLY : init_distrib_icosa_lmdz, transfer_icosa_to_lmdz
49! from dynamico
50  USE domain_mod
51  USE dimensions
52  USE mpi_mod
53  USE mpipara
54  USE disvert_mod
55  USE xios_mod
56  USE time_mod , init_time_icosa=> init_time
57  USE transfert_mod
58 
59! from LMDZ
60  USE mod_grid_phy_lmdz, ONLY : unstructured
61  USE mod_phys_lmdz_omp_data, ONLY: klon_omp
62  USE transfert_mod
63  USE physics_distribution_mod, ONLY : init_physics_distribution
64   
65 
66  IMPLICIT NONE
67  INTEGER  :: ind,i,j,ij,pos
68  REAL(rstd),POINTER :: bounds_lon(:,:)
69  REAL(rstd),POINTER :: bounds_lat(:,:)
70 
71  REAL(rstd),ALLOCATABLE :: latfi(:)
72  REAL(rstd),ALLOCATABLE :: lonfi(:)
73  REAL(rstd),ALLOCATABLE :: airefi(:)
74  REAL(rstd),ALLOCATABLE :: bounds_latfi(:,:)
75  REAL(rstd),ALLOCATABLE :: bounds_lonfi(:,:)
76  REAL(rstd) :: pseudoalt(llm)
77
78  INTEGER :: nbp_phys, nbp_phys_glo
79 
80!$OMP PARALLEL
81    CALL allocate_field(f_bounds_lon,field_t,type_real,6)
82    CALL allocate_field(f_bounds_lat,field_t,type_real,6)
83    CALL allocate_field(f_p,field_t,type_real,llm+1,name="p_in")
84    CALL allocate_field(f_pks,field_t,type_real)
85    CALL allocate_field(f_pk,field_t,type_real,llm)
86    CALL allocate_field(f_p_layer,field_t,type_real,llm,name="p_layer_in")
87    CALL allocate_field(f_theta,field_t,type_real,llm)
88    CALL allocate_field(f_phi,field_t,type_real,llm,name="phi_in")
89    CALL allocate_field(f_Temp,field_t,type_real,llm,name="Temp_in")
90    CALL allocate_field(f_ulon,field_t,type_real,llm,name="ulon_in")
91    CALL allocate_field(f_ulat,field_t,type_real,llm,name="ulat_in")
92    CALL allocate_field(f_dulon,field_t,type_real,llm,name="dulon_out")
93    CALL allocate_field(f_dulat,field_t,type_real,llm,name="dulat_out")
94    CALL allocate_field(f_dTemp,field_t,type_real,llm,name="dTemp_out")
95    CALL allocate_field(f_dq,field_t,type_real,llm,nqtot,name="dq_out")
96    CALL allocate_field(f_dps,field_t,type_real,name="dps_out")
97    CALL allocate_field(f_duc,field_t,type_real,3,llm)   
98
99    CALL init_message(f_dps,req_i0,req_dps0)
100    CALL init_message(f_dulon,req_i0,req_dulon0)
101    CALL init_message(f_dulat,req_i0,req_dulat0)
102    CALL init_message(f_dTemp,req_i0,req_dTemp0)
103    CALL init_message(f_dq,req_i0,req_dq0)
104!$OMP END PARALLEL   
105
106    nbp_phys=0
107    DO ind=1,ndomain
108      CALL swap_dimensions(ind)
109      DO j=jj_begin,jj_end
110        DO i=ii_begin,ii_end
111          IF (domain(ind)%own(i,j)) nbp_phys=nbp_phys+1
112        ENDDO
113      ENDDO
114    ENDDO
115   
116
117!initialize LMDZ5 physic mpi decomposition
118    CALL MPI_ALLREDUCE(nbp_phys,nbp_phys_glo,1,MPI_INTEGER,MPI_SUM,comm_icosa,ierr)
119    CALL init_physics_distribution(unstructured, 6, nbp_phys, 1, nbp_phys_glo, llm, comm_icosa)
120   
121    DO ind=1,ndomain
122        CALL swap_dimensions(ind)
123        CALL swap_geometry(ind)
124        bounds_lon=f_bounds_lon(ind)
125        bounds_lat=f_bounds_lat(ind)
126        DO j=jj_begin,jj_end
127          DO i=ii_begin,ii_end
128            ij=(j-1)*iim+i
129            CALL xyz2lonlat(xyz_v(ij+z_rup,:), bounds_lon(ij,1), bounds_lat(ij,1))
130            CALL xyz2lonlat(xyz_v(ij+z_up,:), bounds_lon(ij,2), bounds_lat(ij,2))
131            CALL xyz2lonlat(xyz_v(ij+z_lup,:), bounds_lon(ij,3), bounds_lat(ij,3))
132            CALL xyz2lonlat(xyz_v(ij+z_ldown,:), bounds_lon(ij,4), bounds_lat(ij,4))
133            CALL xyz2lonlat(xyz_v(ij+z_down,:), bounds_lon(ij,5), bounds_lat(ij,5))
134            CALL xyz2lonlat(xyz_v(ij+z_rdown,:), bounds_lon(ij,6), bounds_lat(ij,6))
135         ENDDO
136       ENDDO           
137    ENDDO
138         
139!$OMP PARALLEL
140    CALL initialize_physics_omp
141!$OMP END PARALLEL           
142
143    CALL xios_set_context   
144
145
146     
147
148  END SUBROUTINE initialize_physics
149
150
151  SUBROUTINE initialize_physics_omp
152  USE distrib_icosa_lmdz_mod, ONLY : init_distrib_icosa_lmdz, transfer_icosa_to_lmdz
153! from dynamico
154  USE domain_mod
155  USE dimensions
156  USE mpi_mod
157  USE mpipara
158  USE disvert_mod
159  USE xios_mod
160  USE time_mod , ONLY: init_time_icosa=> init_time, dt, itaumax, itau_physics
161  USE omp_para
162
163! from LMDZ
164  USE mod_grid_phy_lmdz, ONLY : unstructured
165  USE mod_phys_lmdz_omp_data, ONLY: klon_omp
166  USE time_phylmdz_mod, ONLY: init_time_lmdz => init_time
167  USE transfert_mod
168  USE physics_distribution_mod, ONLY : init_physics_distribution
169  USE dimphy, ONLY: init_dimphy
170  USE geometry_mod, ONLY : init_geometry
171  USE vertical_layers_mod, ONLY : init_vertical_layers
172!  USE planete_mod, ONLY: ini_planete_mod
173  USE cpdet_phy_mod, ONLY: init_cpdet_phy
174  USE infotrac_phy, ONLY: init_infotrac_phy
175 
176  USE netcdf 
177 
178  IMPLICIT NONE
179
180
181
182  INTEGER  :: ind,i,j,k,ij,pos
183  REAL(rstd),POINTER :: bounds_lon(:,:)
184  REAL(rstd),POINTER :: bounds_lat(:,:)
185 
186  REAL(rstd),ALLOCATABLE :: latfi(:)
187  REAL(rstd),ALLOCATABLE :: lonfi(:)
188  REAL(rstd),ALLOCATABLE :: airefi(:)
189  REAL(rstd),ALLOCATABLE :: bounds_latfi(:,:)
190  REAL(rstd),ALLOCATABLE :: bounds_lonfi(:,:)
191  REAL(rstd),ALLOCATABLE :: ind_cell_glo(:)
192  REAL(rstd),ALLOCATABLE :: dx(:)
193  REAL(rstd),ALLOCATABLE :: dy(:)
194
195  REAL(rstd) :: pseudoalt(llm)
196  REAL(rstd) :: aps(llm)
197  REAL(rstd) :: bps(llm)
198  real(rstd) :: scaleheight
199
200  ! For Cp(T)
201  LOGICAL :: cpofT
202  REAL(rstd) :: nu_venus
203  REAL(rstd) :: t0_venus
204
205  ! Calendar related stuff
206  REAL(rstd) :: ptimestep
207  REAL(rstd) :: run_length 
208  INTEGER :: annee_ref 
209  INTEGER :: day_ref   
210  INTEGER :: day_ini
211  INTEGER :: day_end
212  INTEGER :: raz_date
213
214  ! Tracer related information
215  INTEGER :: iflag_trac
216  CHARACTER(len=4)              :: type_trac
217!  CHARACTER(len=30),ALLOCATABLE :: tname(:)    ! tracer short name for restart and diagnostics
218!  CHARACTER(len=33),ALLOCATABLE :: ttext(:)     ! tracer long name for diagnostics
219  TYPE(t_field),POINTER,SAVE    :: f_ind_cell_glo(:)
220 
221  INTEGER :: iflag_phys   
222  INTEGER :: nq
223
224  !! to get values from startphy.nc controle array
225  !! ---------------------------------------------
226  logical :: startphy_file
227  ! NetCDF stuff
228  integer :: status ! NetCDF return code
229  integer :: ncid ! NetCDF file ID
230  integer :: varid ! NetCDF variable ID
231  real :: tab_cntrl(100)
232
233    CALL init_distrib_icosa_lmdz
234   
235    ALLOCATE(latfi(klon_omp))
236    ALLOCATE(lonfi(klon_omp))
237    ALLOCATE(airefi(klon_omp))
238    ALLOCATE(bounds_latfi(klon_omp,6))
239    ALLOCATE(bounds_lonfi(klon_omp,6))
240    ALLOCATE(ind_cell_glo(klon_omp))
241    ALLOCATE(dx(klon_omp))
242    ALLOCATE(dy(klon_omp))
243
244    CALL transfer_icosa_to_lmdz(geom%lat_i,latfi)
245    CALL transfer_icosa_to_lmdz(geom%lon_i,lonfi)
246    CALL transfer_icosa_to_lmdz(f_bounds_lat,bounds_latfi)
247    CALL transfer_icosa_to_lmdz(f_bounds_lon,bounds_lonfi)
248    CALL transfer_icosa_to_lmdz(geom%Ai,airefi)
249
250    CALL allocate_field(f_ind_cell_glo,field_t,type_real)
251   
252    DO ind=1,ndomain
253      IF (.NOT. assigned_domain(ind)  .OR. .NOT. is_omp_level_master ) CYCLE
254      CALL swap_dimensions(ind)
255      CALL swap_geometry(ind)
256      DO j=jj_begin,jj_end
257        DO i=ii_begin,ii_end
258          ij=(j-1)*iim+i
259          f_ind_cell_glo(ind)%rval2d(ij)=domain(ind)%assign_cell_glo(i,j)
260        ENDDO
261      ENDDO
262    ENDDO
263
264     
265    CALL transfer_icosa_to_lmdz(f_ind_cell_glo,ind_cell_glo)
266    CALL deallocate_field(f_ind_cell_glo)
267     
268             
269    ! Initialize dimphy module
270    CALL init_dimphy(klon_omp,llm)
271
272    ! Dummy initializations for dx(),dy() In principle these are not used
273    ! in the physics; but this should be checked further...
274    dx(:)=1 ; dy(:)=1
275    CALL init_geometry(klon_omp,lonfi,latfi,bounds_lonfi,bounds_latfi,&
276                       airefi,INT(ind_cell_glo),dx,dy)
277
278    scaleheight=scale_height/1000. ! Atmospheric scale height (km)
279    aps(1:llm)=0.5*(ap(1:llm)+ap(2:llm+1))
280    bps(1:llm)=0.5*(bp(1:llm)+bp(2:llm+1))
281    pseudoalt(:)=-scaleheight*log(presnivs(:)/preff)
282    CALL init_vertical_layers(llm,preff,scaleheight,ap,bp,aps,bps,presnivs,pseudoalt)
283
284    ! Initialize planet_mod (quite redundant wrt vertical_levels...)
285!    CALL ini_planete_mod(llm,preff,ap,bp)
286
287    ! Initialize tracer names, numbers, etc. for physics
288
289! init tracers model for standard lmdz case
290!$OMP MASTER
291    ALLOCATE(tname(nqtot))
292    ALLOCATE(ttext(nqtot))
293!$OMP END MASTER
294!$OMP BARRIER
295
296! read tname() from traceur.def file
297    IF (is_mpi_root) THEN
298!$OMP MASTER
299    OPEN(unit=42,file="traceur.def",form="formatted",status="old",iostat=ierr)
300    IF (ierr==0) THEN
301      READ(42,*) nq ! should be the same as nqtot
302      IF (nq /= nqtot) THEN
303        WRITE(*,*) "Error: number of tracers in tracer.def should match nqtot!"
304        WRITE(*,*) "       will just use nqtot=",nqtot," tracers"
305      ENDIF
306      DO i=1,nqtot
307        READ(42,*) j,k,tname(i)
308        ttext(i)=trim(tname(i))//"VL1"
309      ENDDO
310      CLOSE(42)
311    ENDIF
312!$OMP END MASTER
313!$OMP BARRIER
314    ENDIF ! of (is_mpi_root)
315
316    DO i=1,nqtot
317      CALL bcast(tname(i))
318      CALL bcast(ttext(i))
319    ENDDO
320
321   ! Get/set some constants for the physics
322
323
324    startphy_file=.true.
325    CALL getin('startphy_file',startphy_file)
326
327    ! value in physics daysec=10087066.76s, important for solar radiation,
328    ! for the rest ptime/timeofday is what matters. To well estimate them,
329    ! day_length must be multiple of dt*itau_physics. Error order:1e-5.
330    day_length=10087000   
331    CALL getin('day_length',day_length)
332   
333    run_length=day_length ! default
334    CALL getin('run_length',run_length)
335   
336    raz_date=0 ! default: 0: no change in date
337    CALL getin('raz_date',raz_date)
338
339    IF (startphy_file) THEN
340      ! Read in some information from the startphy.nc file
341
342      IF (is_mpi_root) THEN
343!$OMP MASTER     
344      status=nf90_open('startphy.nc',NF90_NOWRITE,ncid)
345      if (status.ne.nf90_noerr) then
346        write(*,*)"Failed to open startphy.nc"
347        write(*,*)trim(nf90_strerror(status))
348        stop
349      endif
350
351      status=nf90_inq_varid(ncid,"controle",varid)
352      if (status.ne.nf90_noerr) then
353        write(*,*)"Failed to find controle variable"
354        write(*,*)trim(nf90_strerror(status))
355        stop
356      endif
357
358      status=nf90_get_var(ncid,varid,tab_cntrl)
359      ! extract needed variables from tab_cntrl
360      day_ini=tab_cntrl(13)
361      annee_ref=tab_cntrl(14)
362      ! check if tab_cntrl(1), the stored physics time step
363      ! is the same as the current physics time step (within roundoff precision)
364      if (abs(((itau_physics*dt)-tab_cntrl(1))/(itau_physics*dt))<=1.e-8) then
365        ! Everything OK
366        ptime=modulo((tab_cntrl(15)*tab_cntrl(1))/day_length,1.0)
367      else ! unless raz_date == 1 , we have a problem
368        if (raz_date==1) then
369          ! we reset date to midnight at lon=0 in the physics
370          ptime=0.0
371          day_ini=0
372        else
373          write(*,*)"Error: physics time step in startphy.nc is different"
374          write(*,*)"       from what is specified via run_icosa.def"
375          write(*,*)"       From run_icosa.def:    ",itau_physics*dt
376          write(*,*)"       From startphy.nc file: ",tab_cntrl(1)
377          write(*,*)"       You must reset date to midnight at lon=0"
378          write(*,*)"       by specifying raz_date=1 in your def file"
379          call abort
380        endif
381      endif
382
383      status=nf90_close(ncid)
384!$OMP END MASTER     
385!$OMP BARRIER
386      ENDIF ! of !IF (is_mpi_root)
387     
388      CALL bcast(day_ini)
389      CALL bcast(annee_ref)
390      CALL bcast(ptime)
391
392    ELSE
393      ! required information that is not in tab_cntrl
394      ! has to be default or read from def files
395      day_ini=0
396      annee_ref=1
397      CALL getin('annee_ref',annee_ref)
398      ptime=0.
399    ENDIF
400
401    day_end=day_ini+nint(run_length/day_length)
402
403    ! Other required values which have to be read from def files
404    day_ref=1
405    CALL getin('day_ref',day_ref)
406    iflag_trac=0
407    CALL getin('iflag_trac',iflag_trac)
408
409
410    ! Initialize some physical constants
411    CALL suphec(cpp)
412
413    ! Initialize cpdet_phy module
414! FLAG CP CONSTANT OR NOT... (in LMDZ.COMMON: cpofT in gcm.def)
415! here, it allows to run the physics with a variable Cp together with dynamics
416! with constant Cp.
417! Conversion Theta/T takes this into account (see below)
418    cpofT=.false. ! default, assume Cp to be constant in physics
419    CALL getin('cpofT',cpofT)
420
421    IF (cpofT) THEN
422      nu_venus=0.35
423      t0_venus=460.
424    ELSE
425      nu_venus=0.
426      t0_venus=0.
427    ENDIF
428
429    CALL init_cpdet_phy(cpp,nu_venus,t0_venus)
430 
431    ! Initialize some "temporal and calendar" related variables
432    ptimestep=itau_physics*dt
433    CALL init_time_lmdz(annee_ref,day_ref,day_ini,day_end,ptimestep)
434 
435    ! Initialize tracers in physics
436    CALL init_infotrac_phy(iflag_trac,nqtot,tname,ttext)
437   
438   
439    ! Initializations of some module variables
440   
441!$OMP MASTER
442    ! initialize pday
443    pday = day_ini
444!$OMP END MASTER
445!$OMP BARRIER
446 
447  END SUBROUTINE  initialize_physics_omp
448 
449 
450
451
452  SUBROUTINE physics
453  USE icosa
454  USE time_mod
455  USE disvert_mod
456  USE transfert_mod
457  USE mpipara
458  USE xios_mod
459  USE wxios
460  USE trace
461  USE distrib_icosa_lmdz_mod, ONLY : transfer_icosa_to_lmdz, transfer_lmdz_to_icosa
462  USE physics_external_mod, ONLY : it, f_phis, f_ps, f_theta_rhodz, f_u, f_wflux, f_q
463  USE write_field_mod
464  USE checksum_mod
465! from LMDZ
466  USE mod_phys_lmdz_omp_data, ONLY: klon_omp
467  USE geometry_mod, ONLY : cell_area
468  USE physiq_mod, ONLY: physiq
469  IMPLICIT NONE
470 
471    REAL(rstd),POINTER :: phis(:)
472    REAL(rstd),POINTER :: ps(:)
473    REAL(rstd),POINTER :: theta_rhodz(:,:,:)
474    REAL(rstd),POINTER :: u(:,:)
475    REAL(rstd),POINTER :: wflux(:,:)
476    REAL(rstd),POINTER :: q(:,:,:)
477    REAL(rstd),POINTER :: p(:,:)
478    REAL(rstd),POINTER :: pks(:)
479    REAL(rstd),POINTER :: pk(:,:)
480    REAL(rstd),POINTER :: p_layer(:,:)
481    REAL(rstd),POINTER :: theta(:,:)
482    REAL(rstd),POINTER :: phi(:,:)
483    REAL(rstd),POINTER :: Temp(:,:)
484    REAL(rstd),POINTER :: ulon(:,:)
485    REAL(rstd),POINTER :: ulat(:,:)
486    REAL(rstd),POINTER :: dulon(:,:)
487    REAL(rstd),POINTER :: dulat(:,:)
488    REAL(rstd),POINTER :: dTemp(:,:)
489    REAL(rstd),POINTER :: dq(:,:,:)
490    REAL(rstd),POINTER :: dps(:)
491    REAL(rstd),POINTER :: duc(:,:,:)
492
493
494    INTEGER :: ind,l
495   
496    REAL(rstd),ALLOCATABLE,SAVE :: ps_phy(:)
497!$OMP THREADPRIVATE(ps_phy)
498    REAL(rstd),ALLOCATABLE,SAVE :: p_phy(:,:)
499!$OMP THREADPRIVATE(p_phy)
500    REAL(rstd),ALLOCATABLE,SAVE :: p_layer_phy(:,:)
501!$OMP THREADPRIVATE(p_layer_phy)
502    REAL(rstd),ALLOCATABLE,SAVE :: pk_phy(:,:)
503!$OMP THREADPRIVATE(pk_phy)
504    REAL(rstd),ALLOCATABLE,SAVE :: Temp_phy(:,:)
505!$OMP THREADPRIVATE(Temp_phy)
506    REAL(rstd),ALLOCATABLE,SAVE :: phis_phy(:)
507!$OMP THREADPRIVATE(phis_phy)
508    REAL(rstd),ALLOCATABLE,SAVE :: phi_phy(:,:)
509!$OMP THREADPRIVATE(phi_phy)
510    REAL(rstd),ALLOCATABLE,SAVE :: ulon_phy(:,:)
511!$OMP THREADPRIVATE(ulon_phy)
512    REAL(rstd),ALLOCATABLE,SAVE :: ulat_phy(:,:)
513!$OMP THREADPRIVATE(ulat_phy)
514    REAL(rstd),ALLOCATABLE,SAVE :: q_phy(:,:,:)
515!$OMP THREADPRIVATE(q_phy)
516    REAL(rstd),ALLOCATABLE,SAVE :: wflux_phy(:,:)
517!$OMP THREADPRIVATE(wflux_phy)
518    REAL(rstd),ALLOCATABLE,SAVE :: dulon_phy(:,:)
519!$OMP THREADPRIVATE(dulon_phy)
520    REAL(rstd),ALLOCATABLE,SAVE :: dulat_phy(:,:)
521!$OMP THREADPRIVATE(dulat_phy)
522    REAL(rstd),ALLOCATABLE,SAVE :: dTemp_phy(:,:)
523!$OMP THREADPRIVATE(dTemp_phy)
524    REAL(rstd),ALLOCATABLE,SAVE :: dq_phy(:,:,:)
525!$OMP THREADPRIVATE(dq_phy)
526    REAL(rstd),ALLOCATABLE,SAVE :: dps_phy(:)
527!$OMP THREADPRIVATE(dps_phy)
528    REAL(rstd)   :: dtphy
529    LOGICAL      :: debut
530    LOGICAL      :: lafin
531    LOGICAL,SAVE :: first=.TRUE.
532!$OMP THREADPRIVATE(first)
533   
534    IF(first) THEN
535      debut=.TRUE.
536    ELSE
537      debut=.FALSE.
538    ENDIF
539
540
541    IF(it-itau0>=itaumax) THEN
542      lafin=.TRUE.
543    ELSE
544      lafin=.FALSE.
545    ENDIF
546
547    IF (first) THEN
548      first=.FALSE.
549      CALL init_message(f_u,req_e1_vect,req_u)
550      ALLOCATE(ps_phy(klon_omp))
551      ALLOCATE(p_phy(klon_omp,llm+1))
552      ALLOCATE(p_layer_phy(klon_omp,llm))
553      ALLOCATE(pk_phy(klon_omp,llm))
554      ALLOCATE(Temp_phy(klon_omp,llm))
555      ALLOCATE(phis_phy(klon_omp))
556      ALLOCATE(phi_phy(klon_omp,llm))
557      ALLOCATE(ulon_phy(klon_omp,llm))
558      ALLOCATE(ulat_phy(klon_omp,llm))
559      ALLOCATE(q_phy(klon_omp,llm,nqtot))
560      ALLOCATE(wflux_phy(klon_omp,llm))
561      ALLOCATE(dulon_phy(klon_omp,llm))
562      ALLOCATE(dulat_phy(klon_omp,llm))
563      ALLOCATE(dTemp_phy(klon_omp,llm))
564      ALLOCATE(dq_phy(klon_omp,llm,nqtot))
565      ALLOCATE(dps_phy(klon_omp))
566!$OMP BARRIER
567    ENDIF
568
569
570!$OMP MASTER       
571!    CALL update_calendar(it)
572!$OMP END MASTER
573!$OMP BARRIER
574    dtphy=itau_physics*dt
575   
576   
577   
578    CALL transfert_message(f_u,req_u)
579   
580    DO ind=1,ndomain
581      CALL swap_dimensions(ind)
582      IF (assigned_domain(ind)) THEN
583        CALL swap_geometry(ind)
584     
585        phis=f_phis(ind)
586        ps=f_ps(ind)
587        theta_rhodz=f_theta_rhodz(ind)
588        u=f_u(ind)
589        q=f_q(ind)
590        wflux=f_wflux(ind)
591        p=f_p(ind)
592        pks=f_pks(ind)
593        pk=f_pk(ind)
594        p_layer=f_p_layer(ind)
595        theta=f_theta(ind)
596        phi=f_phi(ind)
597        Temp=f_Temp(ind)
598        ulon=f_ulon(ind)
599        ulat=f_ulat(ind)
600           
601        CALL grid_icosa_to_physics
602
603      ENDIF
604    ENDDO
605   
606!$OMP BARRIER
607!$OMP MASTER
608    CALL SYSTEM_CLOCK(start_clock)
609!$OMP END MASTER
610    CALL trace_start("physic")
611!    CALL trace_off()
612
613
614!    CALL writeField("p_in",f_p)
615!    CALL writeField("p_layer_in",f_p_layer)
616!    CALL writeField("phi_in",f_phi)
617!    CALL writeField("phis_in",f_phis)
618!    CALL writeField("ulon_in",f_ulon)
619!    CALL writeField("ulat_in",f_ulat)
620!    CALL writeField("Temp_in",f_Temp)
621!    CALL writeField("q_in",f_q)
622!    CALL writeField("wflux_in",f_wflux)
623
624!    CALL checksum(f_p)
625!    CALL checksum(f_p_layer)
626!    CALL checksum(f_phi)
627!    CALL checksum(f_phis)
628!    CALL checksum(f_ulon)
629!    CALL checksum(f_ulat)
630!    CALL checksum(f_Temp)
631!    CALL checksum(f_q)
632!    CALL checksum(f_wflux)
633
634    CALL transfer_icosa_to_lmdz(f_p      , p_phy)
635    CALL transfer_icosa_to_lmdz(f_p_layer, p_layer_phy)
636    CALL transfer_icosa_to_lmdz(f_pk     , pk_phy)
637    CALL transfer_icosa_to_lmdz(f_phi    , phi_phy)
638    CALL transfer_icosa_to_lmdz(f_phis   , phis_phy )
639    CALL transfer_icosa_to_lmdz(f_ulon   , ulon_phy )
640    CALL transfer_icosa_to_lmdz(f_ulat   , ulat_phy)
641    CALL transfer_icosa_to_lmdz(f_Temp   , Temp_phy)
642    CALL transfer_icosa_to_lmdz(f_q      , q_phy)
643    CALL transfer_icosa_to_lmdz(f_wflux  , wflux_phy)
644
645    DO l=1,llm
646      ! Warning: In the physics, vertical flux convention is positive if downwards!
647      wflux_phy(:,l)= - wflux_phy(:,l)*cell_area(:)
648      ! Compute relative geopotential
649      phi_phy(:,l)=phi_phy(:,l)-phis_phy(:)
650    ENDDO
651   
652    CALL wxios_set_context()
653
654    ! Update pday and ptime to send to physics
655!$OMP MASTER
656    ptime=ptime+dtphy/day_length
657    IF (ptime >= 1.) THEN
658      ptime=ptime-1
659      pday=pday+1
660    ENDIF
661!$OMP END MASTER
662!$OMP BARRIER   
663
664    CALL physiq(klon_omp, llm, nqtot, &
665                debut, lafin, &
666                pday, ptime, dtphy, &
667                p_phy, p_layer_phy, pk_phy, &
668                phi_phy, phis_phy, presnivs, &
669                ulon_phy, ulat_phy, Temp_phy, q_phy, &
670                wflux_phy, &
671                dulon_phy, dulat_phy, dTemp_phy, dq_phy, dps_phy)
672   
673    CALL transfer_lmdz_to_icosa(dulon_phy, f_dulon )
674    CALL transfer_lmdz_to_icosa(dulat_phy, f_dulat )
675    CALL transfer_lmdz_to_icosa(dTemp_phy, f_dTemp )
676    CALL transfer_lmdz_to_icosa(dq_phy   , f_dq )
677    CALL transfer_lmdz_to_icosa(dps_phy  , f_dps )
678 
679!    CALL writeField("dulon_out",f_dulon)
680!    CALL writeField("dulat_out",f_dulat)
681!    CALL writeField("dTemp_out",f_dTemp)
682!    CALL writeField("dq_out",f_dq)
683!    CALL writeField("dps_out",f_dps)
684
685!    CALL checksum(f_dulon)
686!    CALL checksum(f_dulat)
687!    CALL checksum(f_dTemp)
688!    CALL checksum(f_dq)
689!    CALL checksum(f_dps)
690   
691    CALL send_message(f_dps,req_dps0)
692    CALL send_message(f_dulon,req_dulon0)
693    CALL send_message(f_dulat,req_dulat0)
694    CALL send_message(f_dTemp,req_dTemp0)
695    CALL send_message(f_dq,req_dq0)
696
697    CALL wait_message(req_dps0)
698    CALL wait_message(req_dulon0)
699    CALL wait_message(req_dulat0)
700    CALL wait_message(req_dTemp0)
701    CALL wait_message(req_dq0)
702
703
704!    CALL trace_on()
705    CALL trace_end("physic")
706!$OMP MASTER
707    CALL SYSTEM_CLOCK(stop_clock)
708    count_clock=count_clock+stop_clock-start_clock
709!$OMP END MASTER
710
711!$OMP BARRIER                       
712
713    DO ind=1,ndomain
714      CALL swap_dimensions(ind)
715      IF (assigned_domain(ind)) THEN
716        CALL swap_geometry(ind)
717
718        theta_rhodz=f_theta_rhodz(ind)
719        u=f_u(ind)
720        q=f_q(ind)
721        ps=f_ps(ind)
722        dulon=f_dulon(ind)
723        dulat=f_dulat(ind)
724        Temp=f_temp(ind)
725        dTemp=f_dTemp(ind)
726        dq=f_dq(ind)
727        dps=f_dps(ind)
728        duc=f_duc(ind)
729        p=f_p(ind)
730        pks=f_pks(ind)
731        pk=f_pk(ind)
732     
733        CALL grid_physics_to_icosa
734      ENDIF
735    ENDDO
736
737!$OMP BARRIER
738    CALL xios_set_context   
739   
740 
741  CONTAINS
742
743    SUBROUTINE grid_icosa_to_physics
744    USE pression_mod
745    USE exner_mod
746    USE theta2theta_rhodz_mod
747    USE geopotential_mod
748    USE omp_para
749    USE cpdet_phy_mod, only: tpot2t
750    IMPLICIT NONE
751   
752    REAL(rstd) :: uc(3)
753    INTEGER :: i,j,ij,l
754    REAL    :: tmpT(1),tmptheta(1),tmppk(1)
755   
756
757! compute pression
758
759      DO    l    = ll_begin,ll_endp1
760        DO j=jj_begin,jj_end
761          DO i=ii_begin,ii_end
762            ij=(j-1)*iim+i
763            p(ij,l) = ap(l) + bp(l) * ps(ij)
764          ENDDO
765        ENDDO
766      ENDDO
767
768!$OMP BARRIER
769
770! compute exner
771       
772       IF (is_omp_first_level) THEN
773         DO j=jj_begin,jj_end
774            DO i=ii_begin,ii_end
775               ij=(j-1)*iim+i
776               pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
777            ENDDO
778         ENDDO
779       ENDIF
780
781       ! 3D : pk
782       DO l = ll_begin,ll_end
783          DO j=jj_begin,jj_end
784             DO i=ii_begin,ii_end
785                ij=(j-1)*iim+i
786                pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa
787             ENDDO
788          ENDDO
789       ENDDO
790
791!$OMP BARRIER
792
793!   compute theta, temperature and pression at layer
794! Conversion Theta/T takes Cp(T) into account
795    DO    l    = ll_begin, ll_end
796      DO j=jj_begin,jj_end
797        DO i=ii_begin,ii_end
798          ij=(j-1)*iim+i
799          theta(ij,l) = theta_rhodz(ij,l,1) / ((p(ij,l)-p(ij,l+1))/g)
800          tmptheta(1) = theta(ij,l)
801          tmppk(1)    = pk(ij,l)
802          call tpot2t(1,tmptheta,tmpT,tmppk)
803          Temp(ij,l)  = tmpT(1)
804          p_layer(ij,l)=preff*(pk(ij,l)/cpp)**(1./kappa)
805        ENDDO
806      ENDDO
807    ENDDO
808
809
810!!! Compute geopotential
811!!  This computation (with teta = cp T / pk !) is identical to
812!!     delta phi = R/RMD T/p delta p         (r=R/RMD=cpp*kappa)
813       
814  ! for first layer
815  IF (is_omp_first_level) THEN
816    DO j=jj_begin,jj_end
817      DO i=ii_begin,ii_end
818        ij=(j-1)*iim+i
819        phi( ij,1 ) = phis( ij ) + ( cpp*Temp(ij,1)/pk(ij,1) ) &
820                                 * ( pks(ij) - pk(ij,1) )
821      ENDDO
822    ENDDO
823  ENDIF
824!!-> implicit flush on phi(:,1)
825         
826  ! for other layers
827  DO l = ll_beginp1, ll_end
828    DO j=jj_begin,jj_end
829      DO i=ii_begin,ii_end
830        ij=(j-1)*iim+i
831        phi(ij,l) =  0.5 *cpp*( Temp(ij,l)/pk(ij,l)+Temp(ij,l-1)/pk(ij,l-1) ) &
832                         * ( pk(ij,l-1) - pk(ij,l) )
833      ENDDO
834    ENDDO
835  ENDDO       
836
837!$OMP BARRIER
838
839
840  IF (is_omp_first_level) THEN
841    DO l = 2, llm
842      DO j=jj_begin,jj_end
843! ---> Bug compilo intel ici en openmp
844! ---> Couper la boucle
845       IF (j==jj_end+1) PRINT*,"this message must not be printed"
846        DO i=ii_begin,ii_end
847          ij=(j-1)*iim+i
848          phi(ij,l) = phi(ij,l)+ phi(ij,l-1)
849        ENDDO
850      ENDDO
851    ENDDO
852! --> IMPLICIT FLUSH on phi --> non
853  ENDIF
854
855! compute wind centered lon lat compound
856    DO l=ll_begin,ll_end
857      DO j=jj_begin,jj_end
858        DO i=ii_begin,ii_end
859          ij=(j-1)*iim+i
860          uc(:)=1/Ai(ij)*                                                                                                &
861                        ( ne(ij,right)*u(ij+u_right,l)*le(ij+u_right)*((xyz_v(ij+z_rdown,:)+xyz_v(ij+z_rup,:))/2-centroid(ij,:))  &
862                         + ne(ij,rup)*u(ij+u_rup,l)*le(ij+u_rup)*((xyz_v(ij+z_rup,:)+xyz_v(ij+z_up,:))/2-centroid(ij,:))          &
863                         + ne(ij,lup)*u(ij+u_lup,l)*le(ij+u_lup)*((xyz_v(ij+z_up,:)+xyz_v(ij+z_lup,:))/2-centroid(ij,:))          &
864                         + ne(ij,left)*u(ij+u_left,l)*le(ij+u_left)*((xyz_v(ij+z_lup,:)+xyz_v(ij+z_ldown,:))/2-centroid(ij,:))    &
865                         + ne(ij,ldown)*u(ij+u_ldown,l)*le(ij+u_ldown)*((xyz_v(ij+z_ldown,:)+xyz_v(ij+z_down,:))/2-centroid(ij,:))&
866                         + ne(ij,rdown)*u(ij+u_rdown,l)*le(ij+u_rdown)*((xyz_v(ij+z_down,:)+xyz_v(ij+z_rdown,:))/2-centroid(ij,:)))
867          ulon(ij,l)=sum(uc(:)*elon_i(ij,:))
868          ulat(ij,l)=sum(uc(:)*elat_i(ij,:))
869        ENDDO
870      ENDDO
871    ENDDO
872
873!$OMP BARRIER
874    END SUBROUTINE grid_icosa_to_physics
875
876
877    SUBROUTINE grid_physics_to_icosa
878    USE theta2theta_rhodz_mod
879    USE omp_para
880    USE cpdet_phy_mod, only: t2tpot
881    IMPLICIT NONE
882      INTEGER :: i,j,ij,l,iq
883      REAL    :: tmptheta(1),tmpT(1),tmppk(1)
884         
885      DO l=ll_begin,ll_end
886        DO j=jj_begin,jj_end
887          DO i=ii_begin,ii_end
888            ij=(j-1)*iim+i
889            duc(ij,:,l)=dulon(ij,l)*elon_i(ij,:)+dulat(ij,l)*elat_i(ij,:)
890          ENDDO
891        ENDDO
892      ENDDO
893
894      DO l=ll_begin,ll_end
895        DO j=jj_begin,jj_end
896          DO i=ii_begin,ii_end
897            ij=(j-1)*iim+i
898            u(ij+u_right,l) = u(ij+u_right,l) + dtphy * sum( 0.5*(duc(ij,:,l) + duc(ij+t_right,:,l))*ep_e(ij+u_right,:) )
899            u(ij+u_lup,l) = u(ij+u_lup,l) + dtphy * sum( 0.5*(duc(ij,:,l) + duc(ij+t_lup,:,l))*ep_e(ij+u_lup,:) )
900            u(ij+u_ldown,l) = u(ij+u_ldown,l) + dtphy*sum( 0.5*(duc(ij,:,l) + duc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:) )
901          ENDDO
902        ENDDO
903      ENDDO         
904
905      DO l=ll_begin,ll_end
906        DO j=jj_begin,jj_end
907          DO i=ii_begin,ii_end
908            ij=(j-1)*iim+i
909            Temp(ij,l)=Temp(ij,l)+ dtphy * dTemp(ij,l)
910          ENDDO
911        ENDDO
912      ENDDO         
913     
914      DO iq=1,nqtot
915        DO l=ll_begin,ll_end
916          DO j=jj_begin,jj_end
917            DO i=ii_begin,ii_end
918              ij=(j-1)*iim+i
919              q(ij,l,iq)=q(ij,l,iq)+ dtphy * dq(ij,l,iq)
920            ENDDO
921          ENDDO
922        ENDDO
923      ENDDO
924
925!$OMP BARRIER
926     
927     IF (is_omp_first_level) THEN
928       DO j=jj_begin,jj_end
929         DO i=ii_begin,ii_end
930           ij=(j-1)*iim+i
931           ps(ij)=ps(ij)+ dtphy * dps(ij)
932          ENDDO
933       ENDDO
934     ENDIF
935
936!     CALL compute_temperature2theta_rhodz(ps,Temp,theta_rhodz,0)
937
938! compute pression
939!$OMP BARRIER
940      DO    l    = ll_begin,ll_endp1
941        DO j=jj_begin,jj_end
942          DO i=ii_begin,ii_end
943            ij=(j-1)*iim+i
944            p(ij,l) = ap(l) + bp(l) * ps(ij)
945          ENDDO
946        ENDDO
947      ENDDO
948
949!$OMP BARRIER
950
951! compute exner
952       
953       IF (is_omp_first_level) THEN
954         DO j=jj_begin,jj_end
955            DO i=ii_begin,ii_end
956               ij=(j-1)*iim+i
957               pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
958            ENDDO
959         ENDDO
960       ENDIF
961
962       ! 3D : pk
963       DO l = ll_begin,ll_end
964          DO j=jj_begin,jj_end
965             DO i=ii_begin,ii_end
966                ij=(j-1)*iim+i
967                pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa
968             ENDDO
969          ENDDO
970       ENDDO
971
972!$OMP BARRIER
973
974!   compute theta, temperature and pression at layer
975! Conversion Theta/T takes Cp(T) into account
976    DO    l    = ll_begin, ll_end
977      DO j=jj_begin,jj_end
978        DO i=ii_begin,ii_end
979          ij=(j-1)*iim+i
980          tmpT(1)  = Temp(ij,l)
981          tmppk(1) = pk(ij,l)
982          call t2tpot(1,tmpT,tmptheta,tmppk)
983          theta_rhodz(ij,l,1) = tmptheta(1) * ((p(ij,l)-p(ij,l+1))/g)
984        ENDDO
985      ENDDO
986    ENDDO
987   
988    END SUBROUTINE grid_physics_to_icosa
989
990
991
992  END SUBROUTINE physics
993
994
995
996
997
998END MODULE interface_icosa_lmdz_mod
Note: See TracBrowser for help on using the repository browser.