Ignore:
Timestamp:
Nov 12, 2018, 1:52:29 PM (6 years ago)
Author:
Laurent Fairhead
Message:

Inclusion of Yann's latest (summer/fall 2018) modifications for
convergence of DYNAMICO/LMDZ physics
YM/LF

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/DYNAMICO-conv/libf/phylmd/readaerosol.F90

    r2841 r3413  
    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.