Changeset 3435 for LMDZ6/trunk/libf/misc

Timestamp:

Jan 22, 2019, 4:21:59 PM (5 years ago)

Author:

Laurent Fairhead

Message:

"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

Location:

LMDZ6/trunk/libf/misc

Files:

• handle_err_m.F90 (modified) (1 diff)
• regr1_step_av_m.F90 (added)
• wxios.F90 (modified) (9 diffs)

LMDZ6/trunk/libf/misc/handle_err_m.F90

@@ -39 +39 @@
           end if
        end if
-       call abort_gcm("NetCDF95 handle_err", "", 1)
+       call abort_physic("NetCDF95 handle_err", "", 1)
     end if
 

LMDZ6/trunk/libf/misc/wxios.F90

@@ -15 +15 @@
     
     INTEGER, SAVE :: g_comm
-    CHARACTER(len=100), SAVE :: g_ctx_name
+    CHARACTER(len=100), SAVE :: g_ctx_name ="LMDZ"
     TYPE(xios_context), SAVE :: g_ctx
 !$OMP THREADPRIVATE(g_comm,g_cts_name,g_ctx)
@@ -136 +136 @@
         g_ctx_name = xios_ctx_name
         
-        ! Si couple alors init fait dans cpl_init
-        IF (.not. PRESENT(type_ocean)) THEN
-            CALL wxios_context_init()
-        ENDIF
+!        ! Si couple alors init fait dans cpl_init
+!        IF (.not. PRESENT(type_ocean)) THEN
+!            CALL wxios_context_init()
+!        ENDIF
 
     END SUBROUTINE wxios_init
@@ -145 +145 @@
     SUBROUTINE wxios_context_init()
         USE print_control_mod, ONLY : prt_level, lunout
-!        USE mod_phys_lmdz_mpi_data, ONLY : COMM_LMDZ_PHY
+        USE mod_phys_lmdz_mpi_data, ONLY : COMM_LMDZ_PHY
         IMPLICIT NONE
 
@@ -152 +152 @@
 !$OMP MASTER
         !Initialisation du contexte:
-        CALL xios_context_initialize(g_ctx_name, g_comm)
+        !!CALL xios_context_initialize(g_ctx_name, g_comm)
+        CALL xios_context_initialize(g_ctx_name, COMM_LMDZ_PHY)
         CALL xios_get_handle(g_ctx_name, xios_ctx)    !Récupération
         CALL xios_set_current_context(xios_ctx)            !Activation
@@ -165 +166 @@
 !$OMP END MASTER
     END SUBROUTINE wxios_context_init
+
+
+    SUBROUTINE wxios_set_context()
+        IMPLICIT NONE
+        TYPE(xios_context) :: xios_ctx
+
+       !$OMP MASTER
+        CALL xios_get_handle(g_ctx_name, xios_ctx)    !Récupération
+        CALL xios_set_current_context(xios_ctx)            !Activation
+       !$OMP END MASTER
+
+    END SUBROUTINE wxios_set_context
 
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -203 +216 @@
             CASE DEFAULT
                 abort_message = 'wxios_set_cal: Mauvais choix de calendrier'
-                CALL abort_gcm('Gcm:Xios',abort_message,1)
+                CALL abort_physic('Gcm:Xios',abort_message,1)
         END SELECT
         
@@ -237 +250 @@
     ! Pour initialiser un domaine !!!!!!!!!!!!!!!!!!!!
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    SUBROUTINE wxios_domain_param(dom_id, is_sequential, ni, nj, ni_glo, nj_glo,        &
-                                    ibegin, iend, ii_begin, ii_end, jbegin, jend,       &
-                                    data_ni, data_ibegin, data_iend,                    &
-                                    io_lat, io_lon,is_south_pole,mpi_rank)
-         
-
-        USE print_control_mod, ONLY : prt_level, lunout
-        IMPLICIT NONE
-
+    SUBROUTINE wxios_domain_param(dom_id)
+       USE dimphy, only: klon
+       USE mod_phys_lmdz_transfert_para, ONLY: gather, bcast
+       USE mod_phys_lmdz_para, only: jj_nb, jj_begin, jj_end, ii_begin, ii_end, &
+                                     mpi_size, mpi_rank, klon_mpi, &
+                                     is_sequential, is_south_pole_dyn
+       USE mod_grid_phy_lmdz, only: nbp_lon, nbp_lat, klon_glo         
+       USE print_control_mod, ONLY : prt_level, lunout
+       USE geometry_mod
+
+       IMPLICIT NONE
         CHARACTER(len=*),INTENT(IN) :: dom_id ! domain identifier
-        LOGICAL,INTENT(IN) :: is_sequential ! flag
-        INTEGER,INTENT(IN) :: ni ! local MPI domain number of longitudes
-        INTEGER,INTENT(IN) :: nj ! local MPI domain number of latitudes
-        INTEGER,INTENT(IN) :: ni_glo ! global grid number of longitudes
-        INTEGER,INTENT(IN) :: nj_glo ! global grid number of latitudes
-        INTEGER,INTENT(IN) :: ibegin ! start index, on global grid, of local MPI domain 
-        INTEGER,INTENT(IN) :: iend ! end index, on global grid, of local MPI domain
-        INTEGER,INTENT(IN) :: ii_begin ! i index at which local data starts (first row)
-        INTEGER,INTENT(IN) :: ii_end ! i index at which local data ends (last row)
-        INTEGER,INTENT(IN) :: jbegin ! start index, on global grid, of local MPI domain 
-        INTEGER,INTENT(IN) :: jend ! end index, on global grid, of local MPI domain
-        INTEGER,INTENT(IN) :: data_ni
-        INTEGER,INTENT(IN) :: data_ibegin
-        INTEGER,INTENT(IN) :: data_iend
-        REAL,INTENT(IN) :: io_lat(:) ! latitudes (of global grid)
-        REAL,INTENT(IN) :: io_lon(:) ! longitudes (of global grid)
-        logical,intent(in) :: is_south_pole ! does this process include the south pole?
-        integer,intent(in) :: mpi_rank ! rank of process
-        
+
+        REAL   :: rlat_glo(klon_glo)
+        REAL   :: rlon_glo(klon_glo)
+        REAL   :: io_lat(nbp_lat)
+        REAL   :: io_lon(nbp_lon)
+        LOGICAL :: mask(nbp_lon,jj_nb) !Masque pour les problèmes de recouvrement MPI
         TYPE(xios_domain) :: dom
+        INTEGER :: i
         LOGICAL :: boool
         
-        !Masque pour les problèmes de recouvrement MPI:
-        LOGICAL :: mask(ni,nj)
+
+
+        CALL gather(latitude_deg,rlat_glo) 
+        CALL bcast(rlat_glo)
+        CALL gather(longitude_deg,rlon_glo)
+        CALL bcast(rlon_glo)
+    
+  !$OMP MASTER  
+        io_lat(1)=rlat_glo(1)
+        io_lat(nbp_lat)=rlat_glo(klon_glo)
+        IF ((nbp_lon*nbp_lat) > 1) then
+          DO i=2,nbp_lat-1
+            io_lat(i)=rlat_glo(2+(i-2)*nbp_lon)
+          ENDDO
+        ENDIF
+
+        IF (klon_glo == 1) THEN
+          io_lon(1)=rlon_glo(1)
+        ELSE 
+          io_lon(1:nbp_lon)=rlon_glo(2:nbp_lon+1)
+        ENDIF
+
         
         !On récupère le handle:
         CALL xios_get_domain_handle(dom_id, dom)
         
-        IF (prt_level >= 10) THEN
-          WRITE(lunout,*) "wxios_domain_param: mpirank=",mpi_rank," ni:",ni," ni_glo:", ni_glo, " nj:", nj, " nj_glo:", nj_glo
-          WRITE(lunout,*) "wxios_domain_param: mpirank=",mpi_rank," ibegin:",ibegin," iend:", iend, " jbegin:", jbegin, " jend:", jend
-          WRITE(lunout,*) "wxios_domain_param: mpirank=",mpi_rank," ii_begin:",ii_begin," ii_end:", ii_end
-          WRITE(lunout,*) "wxios_domain_param: mpirank=",mpi_rank," Size io_lon:", SIZE(io_lon(ibegin:iend)), " io_lat:", SIZE(io_lat(jbegin:jend))
-        ENDIF
-        
         !On parametrise le domaine:
-        CALL xios_set_domain_attr_hdl(dom, ni_glo=ni_glo, ibegin=ibegin-1, ni=ni, type="rectilinear")
-        CALL xios_set_domain_attr_hdl(dom, nj_glo=nj_glo, jbegin=jbegin-1, nj=nj, data_dim=2)
-        CALL xios_set_domain_attr_hdl(dom, lonvalue_1d=io_lon(ibegin:iend), latvalue_1d=io_lat(jbegin:jend))
+        CALL xios_set_domain_attr_hdl(dom, ni_glo=nbp_lon, ibegin=0, ni=nbp_lon, type="rectilinear")
+        CALL xios_set_domain_attr_hdl(dom, nj_glo=nbp_lat, jbegin=jj_begin-1, nj=jj_nb, data_dim=2)
+        CALL xios_set_domain_attr_hdl(dom, lonvalue_1d=io_lon(1:nbp_lon), latvalue_1d=io_lat(jj_begin:jj_end))
+        CALL xios_set_domain_attr("dom_out", domain_ref=dom_id)
+
         !On definit un axe de latitudes pour les moyennes zonales
         IF (xios_is_valid_axis("axis_lat")) THEN
-           CALL xios_set_axis_attr( "axis_lat", n_glo=nj_glo, n=nj, begin=jbegin-1, value=io_lat(jbegin:jend))
+           CALL xios_set_axis_attr( "axis_lat", n_glo=nbp_lat, n=jj_nb, begin=jj_begin-1, value=io_lat(jj_begin:jj_end))
         ENDIF
 
@@ -294 +312 @@
             mask(:,:)=.TRUE.
             if (ii_begin>1) mask(1:ii_begin-1,1) = .FALSE.
-            if (ii_end<ni) mask(ii_end+1:ni,nj) = .FALSE.
+            if (ii_end<nbp_lon) mask(ii_end+1:nbp_lon,jj_nb) = .FALSE.
             ! special case for south pole
-            if ((ii_end.eq.1).and.(is_south_pole)) mask(1:ni,nj)=.true.
+            if ((ii_end==1).and.(is_south_pole_dyn)) mask(1:nbp_lon,jj_nb)=.true.
             IF (prt_level >= 10) THEN
               WRITE(lunout,*) "wxios_domain_param: mpirank=",mpi_rank," mask(:,1)=",mask(:,1)
-              WRITE(lunout,*) "wxios_domain_param: mpirank=",mpi_rank," mask(:,nj)=",mask(:,nj)
+              WRITE(lunout,*) "wxios_domain_param: mpirank=",mpi_rank," mask(:,jj_nb)=",mask(:,jj_nb)
             ENDIF
             CALL xios_set_domain_attr_hdl(dom, mask_2d=mask)
@@ -311 +329 @@
             IF (prt_level >= 10) WRITE(lunout,*) "wxios_domain_param: Invalid domain: ", trim(dom_id)
         END IF
+!$OMP END MASTER
+        
     END SUBROUTINE wxios_domain_param
     
+
+    SUBROUTINE wxios_domain_param_unstructured(dom_id)
+        USE geometry_mod, ONLY : longitude, latitude, boundslon, boundslat,ind_cell_glo
+        USE mod_grid_phy_lmdz, ONLY : nvertex, klon_glo
+        USE mod_phys_lmdz_para 
+        USE nrtype, ONLY : PI
+        USE ioipsl_getin_p_mod, ONLY : getin_p
+        IMPLICIT NONE
+        CHARACTER(len=*),INTENT(IN) :: dom_id ! domain identifier
+        REAL :: lon_mpi(klon_mpi)
+        REAL :: lat_mpi(klon_mpi)
+        REAL :: boundslon_mpi(klon_mpi,nvertex)
+        REAL :: boundslat_mpi(klon_mpi,nvertex)
+        INTEGER :: ind_cell_glo_mpi(klon_mpi)
+        TYPE(xios_domain) :: dom
+        LOGICAL :: remap_output
+
+        CALL gather_omp(longitude*180/PI,lon_mpi)
+        CALL gather_omp(latitude*180/PI,lat_mpi)
+        CALL gather_omp(boundslon*180/PI,boundslon_mpi)
+        CALL gather_omp(boundslat*180/PI,boundslat_mpi)
+        CALL gather_omp(ind_cell_glo,ind_cell_glo_mpi)
+        
+        remap_output=.TRUE.
+        CALL getin_p("remap_output",remap_output)
+
+!$OMP MASTER
+        CALL xios_get_domain_handle(dom_id, dom)
+        
+        !On parametrise le domaine:
+        CALL xios_set_attr(dom, ni_glo=klon_glo, ibegin=ij_begin-1, ni=ij_nb, type="unstructured")
+        CALL xios_set_attr(dom, nvertex=nvertex, lonvalue_1d=lon_mpi, latvalue_1d=lat_mpi, &
+                           bounds_lon_1d=TRANSPOSE(boundslon_mpi), bounds_lat_1d=TRANSPOSE(boundslat_mpi) )
+        CALL xios_set_attr(dom, i_index=ind_cell_glo_mpi(:)-1)
+        IF (remap_output) THEN
+          CALL xios_set_domain_attr("dom_out", domain_ref="dom_regular")
+          CALL xios_set_fieldgroup_attr("dom_out", domain_ref="dom_regular")
+          CALL xios_set_fieldgroup_attr("remap_expr", expr="@this_ref")
+        ENDIF
+!$OMP END MASTER
+
+    END SUBROUTINE wxios_domain_param_unstructured
+
+
+
+
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     ! Pour déclarer un axe vertical !!!!!!!!!!!!!!!