source: LMDZ6/trunk/libf/dynphy_lonlat/inigeomphy_mod.F90 @ 3603

Last change on this file since 3603 was 3435, checked in by Laurent Fairhead, 6 years ago

"Historic" :-) commit merging the physics branch used for DYNAMICO with the LMDZ trunk.
The same physics branch can now be used seamlessly with the traditional lon-lat LMDZ
dynamical core and DYNAMICO.
Testing consisted in running a lon-lat LMDZ bucket simulation with the NPv6.1 physics package
with the original trunk sources and the merged sources. Tests were succesful in the sense that
numeric continuity was preserved in the restart files from both simulation. Further tests
included running both versions of the physics codes for one year in a LMDZOR setting in which
the restart files also came out identical.

Caution:

  • as the physics package now manages unstructured grids, grid information needs to be transmitted

to the surface scheme ORCHIDEE. This means that the interface defined in surf_land_orchidee_mod.F90
is only compatible with ORCHIDEE version orchidee2.1 and later versions. If previous versions of
ORCHIDEE need to be used, the CPP key ORCHIDEE_NOUNSTRUCT needs to be set at compilation time.
This is done automatically if makelmdz/makelmdz_fcm are called with the veget orchidee2.0 switch

  • due to a limitation in XIOS, the time at which limit conditions will be read in by DYNAMICO will be

delayed by one physic timestep with respect to the time it is read in by the lon-lat model. This is caused
by the line

IF (MOD(itime-1, lmt_pas) == 0 .OR. (jour_lu /= jour .AND. grid_type /= unstructured)) THEN ! time to read

in limit_read_mod.F90

Work still needed on COSP integration and XML files for DYNAMICO

EM, YM, LF

File size: 9.2 KB
Line 
1!
2! $Id: $
3!
4MODULE inigeomphy_mod
5
6CONTAINS
7
8SUBROUTINE inigeomphy(iim,jjm,nlayer, &
9                     nbp, communicator, &
10                     rlatu,rlatv,rlonu,rlonv,aire,cu,cv)
11  USE mod_grid_phy_lmdz, ONLY: klon_glo,  & ! number of atmospheric columns (on full grid)
12                               regular_lonlat, &  ! regular longitude-latitude grid type
13                               nbp_lon, nbp_lat, nbp_lev
14  USE mod_phys_lmdz_para, ONLY: klon_omp, & ! number of columns (on local omp grid)
15                                klon_omp_begin, & ! start index of local omp subgrid
16                                klon_omp_end, & ! end index of local omp subgrid
17                                klon_mpi_begin ! start indes of columns (on local mpi grid)
18  USE geometry_mod, ONLY : init_geometry
19  USE physics_distribution_mod, ONLY : init_physics_distribution
20  USE regular_lonlat_mod, ONLY : init_regular_lonlat, &
21                                 east, west, north, south, &
22                                 north_east, north_west, &
23                                 south_west, south_east
24  USE mod_interface_dyn_phys, ONLY :  init_interface_dyn_phys
25  USE nrtype, ONLY: pi
26  USE comvert_mod, ONLY: preff, ap, bp, aps, bps, presnivs, &
27                         scaleheight, pseudoalt
28  USE vertical_layers_mod, ONLY: init_vertical_layers
29  IMPLICIT NONE
30
31  ! =======================================================================
32  ! Initialisation of the physical constants and some positional and
33  ! geometrical arrays for the physics
34  ! =======================================================================
35
36  include "iniprint.h"
37
38  INTEGER, INTENT (IN) :: nlayer ! number of atmospheric layers
39  INTEGER, INTENT (IN) :: iim ! number of atmospheric columns along longitudes
40  INTEGER, INTENT (IN) :: jjm ! number of atompsheric columns along latitudes
41  INTEGER, INTENT(IN) :: nbp ! number of physics columns for this MPI process
42  INTEGER, INTENT(IN) :: communicator ! MPI communicator
43  REAL, INTENT (IN) :: rlatu(jjm+1) ! latitudes of the physics grid
44  REAL, INTENT (IN) :: rlatv(jjm) ! latitude boundaries of the physics grid
45  REAL, INTENT (IN) :: rlonv(iim+1) ! longitudes of the physics grid
46  REAL, INTENT (IN) :: rlonu(iim+1) ! longitude boundaries of the physics grid
47  REAL, INTENT (IN) :: aire(iim+1,jjm+1) ! area of the dynamics grid (m2)
48  REAL, INTENT (IN) :: cu((iim+1)*(jjm+1)) ! cu coeff. (u_covariant = cu * u)
49  REAL, INTENT (IN) :: cv((iim+1)*jjm) ! cv coeff. (v_covariant = cv * v)
50
51  INTEGER :: ibegin, iend, offset
52  INTEGER :: i,j,k
53  CHARACTER (LEN=20) :: modname = 'inigeomphy'
54  CHARACTER (LEN=80) :: abort_message
55  REAL :: total_area_phy, total_area_dyn
56
57  ! boundaries, on global grid
58  REAL,ALLOCATABLE :: boundslon_reg(:,:)
59  REAL,ALLOCATABLE :: boundslat_reg(:,:)
60
61  ! global array, on full physics grid:
62  REAL,ALLOCATABLE :: latfi_glo(:)
63  REAL,ALLOCATABLE :: lonfi_glo(:)
64  REAL,ALLOCATABLE :: cufi_glo(:)
65  REAL,ALLOCATABLE :: cvfi_glo(:)
66  REAL,ALLOCATABLE :: airefi_glo(:)
67  REAL,ALLOCATABLE :: boundslonfi_glo(:,:)
68  REAL,ALLOCATABLE :: boundslatfi_glo(:,:)
69
70  ! local arrays, on given MPI/OpenMP domain:
71  REAL,ALLOCATABLE,SAVE :: latfi(:)
72  REAL,ALLOCATABLE,SAVE :: lonfi(:)
73  REAL,ALLOCATABLE,SAVE :: cufi(:)
74  REAL,ALLOCATABLE,SAVE :: cvfi(:)
75  REAL,ALLOCATABLE,SAVE :: airefi(:)
76  REAL,ALLOCATABLE,SAVE :: boundslonfi(:,:)
77  REAL,ALLOCATABLE,SAVE :: boundslatfi(:,:)
78  INTEGER,ALLOCATABLE,SAVE :: ind_cell_glo_fi(:)
79!$OMP THREADPRIVATE (latfi,lonfi,cufi,cvfi,airefi,boundslonfi,boundslatfi,ind_cell_glo_fi)
80
81  ! Initialize Physics distibution and parameters and interface with dynamics
82  IF (iim*jjm>1) THEN ! general 3D case
83    CALL init_physics_distribution(regular_lonlat,4, &
84                                 nbp,iim,jjm+1,nlayer,communicator)
85  ELSE ! For 1D model
86    CALL init_physics_distribution(regular_lonlat,4, &
87                                 1,1,1,nlayer,communicator)
88  ENDIF
89  CALL init_interface_dyn_phys
90 
91  ! init regular global longitude-latitude grid points and boundaries
92  ALLOCATE(boundslon_reg(iim,2))
93  ALLOCATE(boundslat_reg(jjm+1,2))
94 
95  DO i=1,iim
96   boundslon_reg(i,east)=rlonu(i+1)
97   boundslon_reg(i,west)=rlonu(i)
98  ENDDO
99
100  boundslat_reg(1,north)= PI/2
101  boundslat_reg(1,south)= rlatv(1)
102  DO j=2,jjm
103   boundslat_reg(j,north)=rlatv(j-1)
104   boundslat_reg(j,south)=rlatv(j)
105  ENDDO
106  boundslat_reg(jjm+1,north)= rlatv(jjm)
107  boundslat_reg(jjm+1,south)= -PI/2
108
109  ! Write values in module regular_lonlat_mod
110  CALL init_regular_lonlat(iim,jjm+1, rlonv(1:iim), rlatu, &
111                           boundslon_reg, boundslat_reg)
112
113  ! Generate global arrays on full physics grid
114  ALLOCATE(latfi_glo(klon_glo),lonfi_glo(klon_glo))
115  ALLOCATE(cufi_glo(klon_glo),cvfi_glo(klon_glo))
116  ALLOCATE(airefi_glo(klon_glo))
117  ALLOCATE(boundslonfi_glo(klon_glo,4))
118  ALLOCATE(boundslatfi_glo(klon_glo,4))
119
120  IF (klon_glo>1) THEN ! general case
121    ! North pole
122    latfi_glo(1)=rlatu(1)
123    lonfi_glo(1)=0.
124    cufi_glo(1) = cu(1)
125    cvfi_glo(1) = cv(1)
126    boundslonfi_glo(1,north_east)=0
127    boundslatfi_glo(1,north_east)=PI/2
128    boundslonfi_glo(1,north_west)=2*PI
129    boundslatfi_glo(1,north_west)=PI/2
130    boundslonfi_glo(1,south_west)=2*PI
131    boundslatfi_glo(1,south_west)=rlatv(1)
132    boundslonfi_glo(1,south_east)=0
133    boundslatfi_glo(1,south_east)=rlatv(1)
134    DO j=2,jjm
135      DO i=1,iim
136        k=(j-2)*iim+1+i
137        latfi_glo(k)= rlatu(j)
138        lonfi_glo(k)= rlonv(i)
139        cufi_glo(k) = cu((j-1)*(iim+1)+i)
140        cvfi_glo(k) = cv((j-1)*(iim+1)+i)
141        boundslonfi_glo(k,north_east)=rlonu(i)
142        boundslatfi_glo(k,north_east)=rlatv(j-1)
143        boundslonfi_glo(k,north_west)=rlonu(i+1)
144        boundslatfi_glo(k,north_west)=rlatv(j-1)
145        boundslonfi_glo(k,south_west)=rlonu(i+1)
146        boundslatfi_glo(k,south_west)=rlatv(j)
147        boundslonfi_glo(k,south_east)=rlonu(i)
148        boundslatfi_glo(k,south_east)=rlatv(j)
149      ENDDO
150    ENDDO
151    ! South pole
152    latfi_glo(klon_glo)= rlatu(jjm+1)
153    lonfi_glo(klon_glo)= 0.
154    cufi_glo(klon_glo) = cu((iim+1)*jjm+1)
155    cvfi_glo(klon_glo) = cv((iim+1)*jjm-iim)
156    boundslonfi_glo(klon_glo,north_east)= 0
157    boundslatfi_glo(klon_glo,north_east)= rlatv(jjm)
158    boundslonfi_glo(klon_glo,north_west)= 2*PI
159    boundslatfi_glo(klon_glo,north_west)= rlatv(jjm)
160    boundslonfi_glo(klon_glo,south_west)= 2*PI
161    boundslatfi_glo(klon_glo,south_west)= -PI/2
162    boundslonfi_glo(klon_glo,south_east)= 0
163    boundslatfi_glo(klon_glo,south_east)= -Pi/2
164
165    ! build airefi(), mesh area on physics grid
166    CALL gr_dyn_fi(1,iim+1,jjm+1,klon_glo,aire,airefi_glo)
167    ! Poles are single points on physics grid
168    airefi_glo(1)=sum(aire(1:iim,1))
169    airefi_glo(klon_glo)=sum(aire(1:iim,jjm+1))
170
171    ! Sanity check: do total planet area match between physics and dynamics?
172    total_area_dyn=sum(aire(1:iim,1:jjm+1))
173    total_area_phy=sum(airefi_glo(1:klon_glo))
174    IF (total_area_dyn/=total_area_phy) THEN
175      WRITE (lunout, *) 'inigeomphy: planet total surface discrepancy !!!'
176      WRITE (lunout, *) '     in the dynamics total_area_dyn=', total_area_dyn
177      WRITE (lunout, *) '  but in the physics total_area_phy=', total_area_phy
178      IF (abs(total_area_dyn-total_area_phy)>0.00001*total_area_dyn) THEN
179        ! stop here if the relative difference is more than 0.001%
180        abort_message = 'planet total surface discrepancy'
181        CALL abort_gcm(modname, abort_message, 1)
182      ENDIF
183    ENDIF
184  ELSE ! klon_glo==1, running the 1D model
185    ! just copy over input values
186    latfi_glo(1)=rlatu(1)
187    lonfi_glo(1)=rlonv(1)
188    cufi_glo(1)=cu(1)
189    cvfi_glo(1)=cv(1)
190    airefi_glo(1)=aire(1,1)
191    boundslonfi_glo(1,north_east)=rlonu(1)
192    boundslatfi_glo(1,north_east)=PI/2
193    boundslonfi_glo(1,north_west)=rlonu(2)
194    boundslatfi_glo(1,north_west)=PI/2
195    boundslonfi_glo(1,south_west)=rlonu(2)
196    boundslatfi_glo(1,south_west)=rlatv(1)
197    boundslonfi_glo(1,south_east)=rlonu(1)
198    boundslatfi_glo(1,south_east)=rlatv(1)
199  ENDIF ! of IF (klon_glo>1)
200
201!$OMP PARALLEL
202  ! Now generate local lon/lat/cu/cv/area/bounds arrays
203  ALLOCATE(latfi(klon_omp),lonfi(klon_omp),cufi(klon_omp),cvfi(klon_omp))
204  ALLOCATE(airefi(klon_omp))
205  ALLOCATE(boundslonfi(klon_omp,4))
206  ALLOCATE(boundslatfi(klon_omp,4))
207  ALLOCATE(ind_cell_glo_fi(klon_omp))
208
209
210  offset = klon_mpi_begin - 1
211  airefi(1:klon_omp) = airefi_glo(offset+klon_omp_begin:offset+klon_omp_end)
212  cufi(1:klon_omp) = cufi_glo(offset+klon_omp_begin:offset+klon_omp_end)
213  cvfi(1:klon_omp) = cvfi_glo(offset+klon_omp_begin:offset+klon_omp_end)
214  lonfi(1:klon_omp) = lonfi_glo(offset+klon_omp_begin:offset+klon_omp_end)
215  latfi(1:klon_omp) = latfi_glo(offset+klon_omp_begin:offset+klon_omp_end)
216  boundslonfi(1:klon_omp,:) = boundslonfi_glo(offset+klon_omp_begin:offset+klon_omp_end,:)
217  boundslatfi(1:klon_omp,:) = boundslatfi_glo(offset+klon_omp_begin:offset+klon_omp_end,:)
218  ind_cell_glo_fi(1:klon_omp)=(/ (i,i=offset+klon_omp_begin,offset+klon_omp_end) /)
219
220  ! copy over local grid longitudes and latitudes
221  CALL init_geometry(klon_omp,lonfi,latfi,boundslonfi,boundslatfi, &
222                     airefi,ind_cell_glo_fi,cufi,cvfi)
223
224  ! copy over preff , ap(), bp(), etc
225  CALL init_vertical_layers(nlayer,preff,scaleheight, &
226                            ap,bp,aps,bps,presnivs,pseudoalt)
227
228!$OMP END PARALLEL
229
230
231END SUBROUTINE inigeomphy
232
233END MODULE inigeomphy_mod
234
Note: See TracBrowser for help on using the repository browser.