Ignore:
Timestamp:
Jan 22, 2019, 4:21:59 PM (6 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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/readaerosol.F90

    r2841 r3435  
    44
    55  REAL, SAVE :: not_valid=-333.
     6 
     7  INTEGER, SAVE :: nbp_lon_src
     8!$OMP THREADPRIVATE(nbp_lon_src) 
     9  INTEGER, SAVE :: nbp_lat_src
     10!$OMP THREADPRIVATE(nbp_lat_src) 
     11  REAL, ALLOCATABLE, SAVE    :: psurf_interp(:,:)
     12!$OMP THREADPRIVATE(psurf_interp) 
    613
    714CONTAINS
    … …  
    167174
    168175
     176  SUBROUTINE init_aero_fromfile(flag_aerosol)
     177  USE netcdf
     178  USE mod_phys_lmdz_para
     179  USE mod_grid_phy_lmdz, ONLY: grid_type, unstructured
     180  USE xios
     181  IMPLICIT NONE
     182  INTEGER, INTENT(IN) :: flag_aerosol
     183  REAL,ALLOCATABLE :: lat_src(:)
     184  REAL,ALLOCATABLE :: lon_src(:)
     185  CHARACTER(LEN=*),PARAMETER :: file_aerosol='aerosols.nat.nc'
     186  CHARACTER(LEN=*),PARAMETER :: file_so4='so4.nat.nc'
     187  INTEGER :: klev_src
     188  INTEGER :: ierr ,ncid, dimID, varid
     189  REAL :: null_array(0)
     190
     191  IF (flag_aerosol>0 .AND. grid_type==unstructured) THEN
     192 
     193    IF (is_omp_root) THEN
     194 
     195      IF (is_mpi_root) THEN
     196   
     197        IF (nf90_open(TRIM(file_aerosol), NF90_NOWRITE, ncid) /= NF90_NOERR) THEN
     198          CALL check_err( nf90_open(TRIM(file_so4), NF90_NOWRITE, ncid), "pb open "//trim(file_so4) )
     199        ENDIF
     200
     201        ! Read and test longitudes
     202        CALL check_err( nf90_inq_dimid(ncid, "lon", dimID),"pb inq dim lon")
     203        CALL check_err( nf90_inquire_dimension(ncid, dimID, len = nbp_lon_src),"pb inq dim lon")
     204        CALL check_err( nf90_inq_varid(ncid, 'lon', varid),"pb inq lon" )
     205        ALLOCATE(lon_src(nbp_lon_src))
     206        CALL check_err( nf90_get_var(ncid, varid, lon_src(:)),"pb get lon" )
     207
     208        ! Read and test latitudes
     209        CALL check_err( nf90_inq_dimid(ncid, "lat", dimID),"pb inq dim lat")
     210        CALL check_err( nf90_inquire_dimension(ncid, dimID, len = nbp_lat_src),"pb inq dim lat")
     211        CALL check_err( nf90_inq_varid(ncid, 'lat', varid),"pb inq lat" )
     212        ALLOCATE(lat_src(nbp_lat_src))
     213        CALL check_err( nf90_get_var(ncid, varid, lat_src(:)),"pb get lat" )
     214        IF (nf90_inq_dimid(ncid, 'lev', dimid) /= NF90_NOERR) THEN
     215          IF (nf90_inq_dimid(ncid, 'presnivs', dimid)/= NF90_NOERR) THEN
     216             CALL check_err(nf90_inq_dimid(ncid, 'PRESNIVS', dimid),'dimension lev,PRESNIVS or presnivs not in file')
     217          ENDIF
     218        ENDIF
     219        CALL check_err( nf90_inquire_dimension(ncid, dimid, len = klev_src),"pb inq dim for PRESNIVS or lev" )
     220        CALL check_err( nf90_close(ncid),"pb in close" )   
     221      ENDIF
     222
     223      CALL bcast_mpi(nbp_lat_src)
     224      CALL bcast_mpi(nbp_lon_src)
     225      CALL bcast_mpi(klev_src)
     226
     227      IF (is_mpi_root ) THEN
     228        CALL xios_set_domain_attr("domain_aerosol",nj_glo=nbp_lat_src, nj=nbp_lat_src, jbegin=0, latvalue_1d=lat_src)
     229        CALL xios_set_domain_attr("domain_aerosol",ni_glo=nbp_lon_src, ni=nbp_lon_src, ibegin=0, lonvalue_1d=lon_src)
     230      ELSE
     231        CALL xios_set_domain_attr("domain_aerosol",nj_glo=nbp_lat_src, nj=0, jbegin=0, latvalue_1d=null_array )
     232        CALL xios_set_domain_attr("domain_aerosol",ni_glo=nbp_lon_src, ni=0, ibegin=0, lonvalue_1d=null_array)
     233      ENDIF
     234      CALL xios_set_axis_attr("axis_aerosol",n_glo=klev_src)
     235      CALL xios_set_fieldgroup_attr("aerosols", enabled=.TRUE.)
     236 
     237    ENDIF
     238   
     239  ENDIF   
     240 
     241  END SUBROUTINE init_aero_fromfile
     242
     243
     244   
     245
    169246  SUBROUTINE get_aero_fromfile(varname, cyr, filename, klev_src, pt_ap, pt_b, pt_year, psurf_out, load_out)
    170247!****************************************************************************************
    … …  
    187264    USE netcdf
    188265    USE dimphy
    189     USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, klon_glo, &
    190                                  grid2Dto1D_glo
     266    USE mod_grid_phy_lmdz, ONLY: nbp_lon_=>nbp_lon, nbp_lat_=>nbp_lat, klon_glo, &
     267                                 grid2Dto1D_glo, grid_type, unstructured
    191268    USE mod_phys_lmdz_para
    192269    USE iophy, ONLY : io_lon, io_lat
    193270    USE print_control_mod, ONLY: lunout
     271    USE xios
    194272
    195273    IMPLICIT NONE
    … …  
    205283    REAL, POINTER, DIMENSION(:)           :: pt_b         ! Pointer for describing the vertical levels     
    206284    REAL, POINTER, DIMENSION(:,:,:)       :: pt_year      ! Pointer-variabale from file, 12 month, grid : klon,klev_src
     285    REAL, POINTER, DIMENSION(:,:,:)       :: pt_year_mpi  ! Pointer-variabale from file, 12 month, grid : klon,klev_src
    207286    REAL, DIMENSION(klon,12), INTENT(OUT) :: psurf_out    ! Surface pression for 12 months
     287    REAL, DIMENSION(klon_mpi,12)          :: psurf_out_mpi    ! Surface pression for 12 months
    208288    REAL, DIMENSION(klon,12), INTENT(OUT) :: load_out     ! Aerosol mass load in each column
     289    REAL, DIMENSION(klon_mpi,12)          :: load_out_mpi     ! Aerosol mass load in each column
    209290    INTEGER                               :: nbr_tsteps   ! number of month in file read
    210291
    … …  
    220301    REAL, ALLOCATABLE, DIMENSION(:)       :: varktmp
    221302
    222     REAL, DIMENSION(nbp_lon,nbp_lat,12)   :: psurf_glo2D   ! Surface pression for 12 months on dynamics global grid
     303    REAL,  ALLOCATABLE                    :: psurf_glo2D(:,:,:)   ! Surface pression for 12 months on dynamics global grid
    223304    REAL, DIMENSION(klon_glo,12)          :: psurf_glo1D   ! -"- on physical global grid
    224     REAL, DIMENSION(nbp_lon,nbp_lat,12)   :: load_glo2D    ! Load for 12 months on dynamics global grid
     305    REAL,  ALLOCATABLE                    :: load_glo2D(:,:,:)    ! Load for 12 months on dynamics global grid
    225306    REAL, DIMENSION(klon_glo,12)          :: load_glo1D    ! -"- on physical global grid
    226     REAL, DIMENSION(nbp_lon,nbp_lat)      :: vartmp
    227     REAL, DIMENSION(nbp_lon)              :: lon_src              ! longitudes in file
    228     REAL, DIMENSION(nbp_lat)              :: lat_src, lat_src_inv ! latitudes in file
     307    REAL, ALLOCATABLE, DIMENSION(:,:)     :: vartmp
     308    REAL, ALLOCATABLE,DIMENSION(:)        :: lon_src              ! longitudes in file
     309    REAL, ALLOCATABLE,DIMENSION(:)        :: lat_src, lat_src_inv ! latitudes in file
    229310    LOGICAL                               :: new_file             ! true if new file format detected
    230311    LOGICAL                               :: invert_lat           ! true if the field has to be inverted for latitudes
    231 
    232 
     312    INTEGER                               :: nbp_lon, nbp_lat
     313    LOGICAL,SAVE                          :: first=.TRUE.
     314!$OMP THREADPRIVATE(first)
     315   
     316    IF (grid_type==unstructured) THEN
     317      nbp_lon=nbp_lon_src
     318      nbp_lat=nbp_lat_src
     319    ELSE
     320      nbp_lon=nbp_lon_
     321      nbp_lat=nbp_lat_
     322    ENDIF
     323   
     324    IF (is_mpi_root) THEN
     325   
     326      ALLOCATE(psurf_glo2D(nbp_lon,nbp_lat,12))
     327      ALLOCATE(load_glo2D(nbp_lon,nbp_lat,12))
     328      ALLOCATE(vartmp(nbp_lon,nbp_lat))
     329      ALLOCATE(lon_src(nbp_lon))
     330      ALLOCATE(lat_src(nbp_lat))
     331      ALLOCATE(lat_src_inv(nbp_lat))
     332    ELSE
     333      ALLOCATE(varyear(0,0,0,0))
     334      ALLOCATE(psurf_glo2D(0,0,0))
     335      ALLOCATE(load_glo2D(0,0,0))
     336    ENDIF
     337           
    233338    ! Deallocate pointers
    234339    IF (ASSOCIATED(pt_ap)) DEALLOCATE(pt_ap)
    … …  
    245350       CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, ncid), "pb open "//trim(fname) )
    246351
     352
     353       IF (grid_type/=unstructured) THEN
     354
    247355! Test for equal longitudes and latitudes in file and model
    248356!****************************************************************************************
    249357       ! Read and test longitudes
    250        CALL check_err( nf90_inq_varid(ncid, 'lon', varid),"pb inq lon" )
    251        CALL check_err( nf90_get_var(ncid, varid, lon_src(:)),"pb get lon" )
     358         CALL check_err( nf90_inq_varid(ncid, 'lon', varid),"pb inq lon" )
     359         CALL check_err( nf90_get_var(ncid, varid, lon_src(:)),"pb get lon" )
    252360       
    253        IF (maxval(ABS(lon_src - io_lon)) > 0.001) THEN
    254           WRITE(lunout,*) 'Problem in longitudes read from file : ',TRIM(fname)
    255           WRITE(lunout,*) 'longitudes in file ', TRIM(fname),' : ', lon_src
    256           WRITE(lunout,*) 'longitudes in model :', io_lon
     361         IF (maxval(ABS(lon_src - io_lon)) > 0.001) THEN
     362            WRITE(lunout,*) 'Problem in longitudes read from file : ',TRIM(fname)
     363            WRITE(lunout,*) 'longitudes in file ', TRIM(fname),' : ', lon_src
     364            WRITE(lunout,*) 'longitudes in model :', io_lon
    257365         
    258           CALL abort_physic('get_aero_fromfile', 'longitudes are not the same in file and model',1)
    259        END IF
    260 
    261        ! Read and test latitudes
    262        CALL check_err( nf90_inq_varid(ncid, 'lat', varid),"pb inq lat" )
    263        CALL check_err( nf90_get_var(ncid, varid, lat_src(:)),"pb get lat" )
    264 
    265        ! Invert source latitudes
    266        DO j = 1, nbp_lat
    267           lat_src_inv(j) = lat_src(nbp_lat +1 -j)
    268        END DO
    269 
    270        IF (maxval(ABS(lat_src - io_lat)) < 0.001) THEN
    271           ! Latitudes are the same
    272           invert_lat=.FALSE.
    273        ELSE IF (maxval(ABS(lat_src_inv - io_lat)) < 0.001) THEN
    274           ! Inverted source latitudes correspond to model latitudes
    275           WRITE(lunout,*) 'latitudes will be inverted for file : ',TRIM(fname)
    276           invert_lat=.TRUE.
    277        ELSE
    278           WRITE(lunout,*) 'Problem in latitudes read from file : ',TRIM(fname)
    279           WRITE(lunout,*) 'latitudes in file ', TRIM(fname),' : ', lat_src     
    280           WRITE(lunout,*) 'latitudes in model :', io_lat
    281           CALL abort_physic('get_aero_fromfile', 'latitudes do not correspond between file and model',1)
    282        END IF
    283 
     366            CALL abort_physic('get_aero_fromfile', 'longitudes are not the same in file and model',1)
     367         END IF
     368
     369         ! Read and test latitudes
     370         CALL check_err( nf90_inq_varid(ncid, 'lat', varid),"pb inq lat" )
     371         CALL check_err( nf90_get_var(ncid, varid, lat_src(:)),"pb get lat" )
     372
     373         ! Invert source latitudes
     374         DO j = 1, nbp_lat
     375            lat_src_inv(j) = lat_src(nbp_lat +1 -j)
     376         END DO
     377
     378         IF (maxval(ABS(lat_src - io_lat)) < 0.001) THEN
     379            ! Latitudes are the same
     380            invert_lat=.FALSE.
     381         ELSE IF (maxval(ABS(lat_src_inv - io_lat)) < 0.001) THEN
     382            ! Inverted source latitudes correspond to model latitudes
     383            WRITE(lunout,*) 'latitudes will be inverted for file : ',TRIM(fname)
     384            invert_lat=.TRUE.
     385         ELSE
     386            WRITE(lunout,*) 'Problem in latitudes read from file : ',TRIM(fname)
     387            WRITE(lunout,*) 'latitudes in file ', TRIM(fname),' : ', lat_src     
     388            WRITE(lunout,*) 'latitudes in model :', io_lat
     389            CALL abort_physic('get_aero_fromfile', 'latitudes do not correspond between file and model',1)
     390         END IF
     391       ENDIF
    284392
    285393! 2) Check if old or new file is avalabale.
    … …  
    487595       END IF
    488596
    489 !     - Invert latitudes if necessary
    490        DO imth=1, 12
    491           IF (invert_lat) THEN
    492 
    493              ! Invert latitudes for the variable
    494              varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly
    495              DO k=1,klev_src
    496                 DO j=1,nbp_lat
    497                    DO i=1,nbp_lon
    498                       varyear(i,j,k,imth) = varmth(i,nbp_lat+1-j,k)
    499                    END DO
    500                 END DO
    501              END DO
     597
     598       IF (grid_type/=unstructured) THEN
     599  !     - Invert latitudes if necessary
     600         DO imth=1, 12
     601            IF (invert_lat) THEN
     602
     603               ! Invert latitudes for the variable
     604               varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly
     605               DO k=1,klev_src
     606                  DO j=1,nbp_lat
     607                     DO i=1,nbp_lon
     608                        varyear(i,j,k,imth) = varmth(i,nbp_lat+1-j,k)
     609                     END DO
     610                  END DO
     611               END DO
    502612             
    503              ! Invert latitudes for surface pressure
    504              vartmp(:,:) = psurf_glo2D(:,:,imth)
    505              DO j=1,nbp_lat
    506                 DO i=1,nbp_lon
    507                    psurf_glo2D(i,j,imth)= vartmp(i,nbp_lat+1-j)
    508                 END DO
    509              END DO
     613               ! Invert latitudes for surface pressure
     614               vartmp(:,:) = psurf_glo2D(:,:,imth)
     615               DO j=1,nbp_lat
     616                  DO i=1,nbp_lon
     617                     psurf_glo2D(i,j,imth)= vartmp(i,nbp_lat+1-j)
     618                  END DO
     619               END DO
    510620             
    511              ! Invert latitudes for the load
    512              vartmp(:,:) = load_glo2D(:,:,imth)
    513              DO j=1,nbp_lat
    514                 DO i=1,nbp_lon
    515                    load_glo2D(i,j,imth)= vartmp(i,nbp_lat+1-j)
    516                 END DO
    517              END DO
    518           END IF ! invert_lat
     621               ! Invert latitudes for the load
     622               vartmp(:,:) = load_glo2D(:,:,imth)
     623               DO j=1,nbp_lat
     624                  DO i=1,nbp_lon
     625                     load_glo2D(i,j,imth)= vartmp(i,nbp_lat+1-j)
     626                  END DO
     627               END DO
     628            END IF ! invert_lat
    519629             
    520           ! Do zonal mead at poles and distribut at whole first and last latitude
    521           DO k=1, klev_src
    522              npole=0.  ! North pole, j=1
    523              spole=0.  ! South pole, j=nbp_lat       
    524              DO i=1,nbp_lon
    525                 npole = npole + varyear(i,1,k,imth)
    526                 spole = spole + varyear(i,nbp_lat,k,imth)
    527              END DO
    528              npole = npole/REAL(nbp_lon)
    529              spole = spole/REAL(nbp_lon)
    530              varyear(:,1,    k,imth) = npole
    531              varyear(:,nbp_lat,k,imth) = spole
    532           END DO
    533        END DO ! imth
     630            ! Do zonal mead at poles and distribut at whole first and last latitude
     631            DO k=1, klev_src
     632               npole=0.  ! North pole, j=1
     633               spole=0.  ! South pole, j=nbp_lat       
     634               DO i=1,nbp_lon
     635                  npole = npole + varyear(i,1,k,imth)
     636                  spole = spole + varyear(i,nbp_lat,k,imth)
     637               END DO
     638               npole = npole/REAL(nbp_lon)
     639               spole = spole/REAL(nbp_lon)
     640               varyear(:,1,    k,imth) = npole
     641               varyear(:,nbp_lat,k,imth) = spole
     642            END DO
     643         END DO ! imth
    534644       
    535        ALLOCATE(varyear_glo1D(klon_glo, klev_src, 12), stat=ierr)
    536        IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 3',1)
     645         ALLOCATE(varyear_glo1D(klon_glo, klev_src, 12), stat=ierr)
     646         IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 3',1)
    537647       
    538        ! Transform from 2D to 1D field
    539        CALL grid2Dto1D_glo(varyear,varyear_glo1D)
    540        CALL grid2Dto1D_glo(psurf_glo2D,psurf_glo1D)
    541        CALL grid2Dto1D_glo(load_glo2D,load_glo1D)
    542 
     648         ! Transform from 2D to 1D field
     649         CALL grid2Dto1D_glo(varyear,varyear_glo1D)
     650         CALL grid2Dto1D_glo(psurf_glo2D,psurf_glo1D)
     651         CALL grid2Dto1D_glo(load_glo2D,load_glo1D)
     652     
     653      ENDIF
     654     
    543655    ELSE
    544       ALLOCATE(varyear_glo1D(0,0,0))       
     656        ALLOCATE(varyear_glo1D(0,0,0))       
    545657    END IF ! is_mpi_root .AND. is_omp_root
    546658
    … …  
    566678    IF (ASSOCIATED(pt_year)) DEALLOCATE(pt_year)
    567679    ALLOCATE(pt_year(klon, klev_src, 12), stat=ierr)
     680    ALLOCATE(pt_year_mpi(klon_mpi, klev_src, 12), stat=ierr)
    568681    IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 5',1)
    569682
    570     ! Scatter global field to local domain at local process
    571     CALL scatter(varyear_glo1D, pt_year)
    572     CALL scatter(psurf_glo1D, psurf_out)
    573     CALL scatter(load_glo1D,  load_out)
    574 
     683    IF (grid_type==unstructured) THEN
     684      IF (is_omp_master) THEN
     685        CALL xios_send_field(TRIM(varname)//"_in",varyear)
     686        CALL xios_recv_field(TRIM(varname)//"_out",pt_year_mpi)
     687        CALL xios_send_field("load_"//TRIM(varname)//"_in",load_glo2D)
     688        CALL xios_recv_field("load_"//TRIM(varname)//"_out",load_out_mpi)
     689        IF (first) THEN
     690          ALLOCATE(psurf_interp(klon_mpi,12))
     691          CALL xios_send_field("psurf_aerosol_in",psurf_glo2D)
     692          CALL xios_recv_field("psurf_aerosol_out",psurf_interp)
     693        ENDIF
     694      ENDIF
     695      CALL scatter_omp(pt_year_mpi,pt_year)
     696      CALL scatter_omp(load_out_mpi,load_out)
     697      CALL scatter_omp(psurf_interp,psurf_out)
     698      first=.FALSE.
     699    ELSE
     700      ! Scatter global field to local domain at local process
     701      CALL scatter(varyear_glo1D, pt_year)
     702      CALL scatter(psurf_glo1D, psurf_out)
     703      CALL scatter(load_glo1D,  load_out)
     704    ENDIF
    575705! 7) Test for negative values
    576706!****************************************************************************************
Note: See TracChangeset for help on using the changeset viewer.