Ignore:
Timestamp:
Jul 24, 2024, 4:39:59 PM (11 months ago)
Author:
abarral
Message:

Replace iniprint.h by lmdz_iniprint.f90
(lint) along the way

Location:
LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/calfis.f90

    r5117 r5118  
    3535  USE comconst_mod, ONLY: cpp, daysec, dtphys, dtvr, kappa, pi
    3636  USE comvert_mod, ONLY: preff, presnivs
     37  USE lmdz_iniprint, ONLY: lunout, prt_level
    3738
    3839  IMPLICIT NONE
    … …  
    9798
    9899  include "comgeom2.h"
    99   include "iniprint.h"
    100100
    101101  !    Arguments :

  • LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/inigeomphy_mod.F90

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

  • LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/lmdz_calfis_loc.F90

    r5117 r5118  
    5050    USE comvert_mod, ONLY: preff, presnivs
    5151    USE comconst_mod, ONLY: cpp, daysec, dtphys, dtvr, kappa, pi
     52    USE lmdz_iniprint, ONLY: lunout, prt_level
    5253
    5354    !=======================================================================
    … …  
    112113
    113114    include "comgeom2.h"
    114     include "iniprint.h"
    115115    !    Arguments :
    116116    !    -----------

  • LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/phydev/iniphysiq_mod.F90

    r5117 r5118  
    1 
    21! $Id: iniphysiq.F 1403 2010-07-01 09:02:53Z fairhead $
    32
    … …  
    65CONTAINS
    76
    8 SUBROUTINE iniphysiq(iim,jjm,nlayer, &
    9                      nbp, communicator, &
    10                      punjours, pdayref,ptimestep, &
    11                      rlatu,rlatv,rlonu,rlonv,aire,cu,cv, &
    12                      prad,pg,pr,pcpp,iflag_phys)
    13   USE dimphy, ONLY: init_dimphy
    14   USE inigeomphy_mod, ONLY: inigeomphy
    15   USE lmdz_phys_para, ONLY: klon_omp ! number of columns (on local omp grid)
    16   USE infotrac, ONLY: nqtot, type_trac
    17   USE infotrac_phy, ONLY: init_infotrac_phy
    18   USE inifis_mod, ONLY: inifis
    19   USE phyaqua_mod, ONLY: iniaqua
    20   USE lmdz_physical_constants, ONLY: pi
    21   IMPLICIT NONE
     7  SUBROUTINE iniphysiq(iim, jjm, nlayer, &
     8          nbp, communicator, &
     9          punjours, pdayref, ptimestep, &
     10          rlatu, rlatv, rlonu, rlonv, aire, cu, cv, &
     11          prad, pg, pr, pcpp, iflag_phys)
     12    USE dimphy, ONLY: init_dimphy
     13    USE inigeomphy_mod, ONLY: inigeomphy
     14    USE lmdz_phys_para, ONLY: klon_omp ! number of columns (on local omp grid)
     15    USE infotrac, ONLY: nqtot, type_trac
     16    USE infotrac_phy, ONLY: init_infotrac_phy
     17    USE inifis_mod, ONLY: inifis
     18    USE phyaqua_mod, ONLY: iniaqua
     19    USE lmdz_physical_constants, ONLY: pi
     20    USE lmdz_iniprint, ONLY: lunout, prt_level
     21    IMPLICIT NONE
    2222
    23   !=======================================================================
    24   !   Initialisation of the physical constants and some positional and
    25   !   geometrical arrays for the physics
    26   !=======================================================================
    27  
    28  
    29   include "iniprint.h"
     23    !=======================================================================
     24    !   Initialisation of the physical constants and some positional and
     25    !   geometrical arrays for the physics
     26    !=======================================================================
    3027
    31   REAL,INTENT(IN) :: prad ! radius of the planet (m)
    32   REAL,INTENT(IN) :: pg ! gravitational acceleration (m/s2)
    33   REAL,INTENT(IN) :: pr ! reduced gas constant R/mu
    34   REAL,INTENT(IN) :: pcpp ! specific heat Cp
    35   REAL,INTENT(IN) :: punjours ! length (in s) of a standard day
    36   INTEGER, INTENT (IN) :: nlayer ! number of atmospheric layers
    37   INTEGER, INTENT (IN) :: iim ! number of atmospheric coulumns along longitudes
    38   INTEGER, INTENT (IN) :: jjm  ! number of atompsheric columns along latitudes
    39   INTEGER, INTENT(IN) :: nbp ! number of physics columns for this MPI process
    40   INTEGER, INTENT(IN) :: communicator ! MPI communicator
    41   REAL, INTENT (IN) :: rlatu(jjm+1) ! latitudes of the physics grid
    42   REAL, INTENT (IN) :: rlatv(jjm) ! latitude boundaries of the physics grid
    43   REAL, INTENT (IN) :: rlonv(iim+1) ! longitudes of the physics grid
    44   REAL, INTENT (IN) :: rlonu(iim+1) ! longitude boundaries of the physics grid
    45   REAL, INTENT (IN) :: aire(iim+1,jjm+1) ! area of the dynamics grid (m2)
    46   REAL, INTENT (IN) :: cu((iim+1)*(jjm+1)) ! cu coeff. (u_covariant = cu * u)
    47   REAL, INTENT (IN) :: cv((iim+1)*jjm) ! cv coeff. (v_covariant = cv * v)
    48   INTEGER, INTENT (IN) :: pdayref ! reference day of for the simulation
    49   REAL,INTENT(IN) :: ptimestep !physics time step (s)
    50   INTEGER,INTENT(IN) :: iflag_phys ! type of physics to be called
     28    REAL, INTENT(IN) :: prad ! radius of the planet (m)
     29    REAL, INTENT(IN) :: pg ! gravitational acceleration (m/s2)
     30    REAL, INTENT(IN) :: pr ! reduced gas constant R/mu
     31    REAL, INTENT(IN) :: pcpp ! specific heat Cp
     32    REAL, INTENT(IN) :: punjours ! length (in s) of a standard day
     33    INTEGER, INTENT (IN) :: nlayer ! number of atmospheric layers
     34    INTEGER, INTENT (IN) :: iim ! number of atmospheric coulumns along longitudes
     35    INTEGER, INTENT (IN) :: jjm  ! number of atompsheric columns along latitudes
     36    INTEGER, INTENT(IN) :: nbp ! number of physics columns for this MPI process
     37    INTEGER, INTENT(IN) :: communicator ! MPI communicator
     38    REAL, INTENT (IN) :: rlatu(jjm + 1) ! latitudes of the physics grid
     39    REAL, INTENT (IN) :: rlatv(jjm) ! latitude boundaries of the physics grid
     40    REAL, INTENT (IN) :: rlonv(iim + 1) ! longitudes of the physics grid
     41    REAL, INTENT (IN) :: rlonu(iim + 1) ! longitude boundaries of the physics grid
     42    REAL, INTENT (IN) :: aire(iim + 1, jjm + 1) ! area of the dynamics grid (m2)
     43    REAL, INTENT (IN) :: cu((iim + 1) * (jjm + 1)) ! cu coeff. (u_covariant = cu * u)
     44    REAL, INTENT (IN) :: cv((iim + 1) * jjm) ! cv coeff. (v_covariant = cv * v)
     45    INTEGER, INTENT (IN) :: pdayref ! reference day of for the simulation
     46    REAL, INTENT(IN) :: ptimestep !physics time step (s)
     47    INTEGER, INTENT(IN) :: iflag_phys ! type of physics to be called
    5148
    52   INTEGER :: ibegin,iend,offset
    53   INTEGER :: i,j,k
    54   CHARACTER (LEN=20) :: modname='iniphysiq'
    55   CHARACTER (LEN=80) :: abort_message
     49    INTEGER :: ibegin, iend, offset
     50    INTEGER :: i, j, k
     51    CHARACTER (LEN = 20) :: modname = 'iniphysiq'
     52    CHARACTER (LEN = 80) :: abort_message
    5653
    5754
    58   ! --> initialize physics distribution, global fields and geometry
    59   ! (i.e. things in phy_common or dynphy_lonlat)
    60   CALL inigeomphy(iim,jjm,nlayer, &
    61                nbp, communicator, &
    62                rlatu,rlatv, &
    63                rlonu,rlonv, &
    64                aire,cu,cv)
     55    ! --> initialize physics distribution, global fields and geometry
     56    ! (i.e. things in phy_common or dynphy_lonlat)
     57    CALL inigeomphy(iim, jjm, nlayer, &
     58            nbp, communicator, &
     59            rlatu, rlatv, &
     60            rlonu, rlonv, &
     61            aire, cu, cv)
    6562
    66   ! --> now initialize things specific to the phydev physics package
     63    ! --> now initialize things specific to the phydev physics package
    6764
    68 !$OMP PARALLEL
     65    !$OMP PARALLEL
    6966
    70   ! Initialize physical constants in physics:
    71   CALL inifis(prad,pg,pr,pcpp)
    72  
    73   ! Initialize tracer names, numbers, etc. for physics
    74   CALL init_infotrac_phy(nqtot,type_trac)
     67    ! Initialize physical constants in physics:
     68    CALL inifis(prad, pg, pr, pcpp)
    7569
    76   ! Additional initializations for aquaplanets
    77   IF (iflag_phys>=100) THEN
    78     CALL iniaqua(klon_omp,iflag_phys)
    79   ENDIF
    80 !$OMP END PARALLEL
     70    ! Initialize tracer names, numbers, etc. for physics
     71    CALL init_infotrac_phy(nqtot, type_trac)
    8172
    82 END SUBROUTINE iniphysiq
     73    ! Additional initializations for aquaplanets
     74    IF (iflag_phys>=100) THEN
     75      CALL iniaqua(klon_omp, iflag_phys)
     76    ENDIF
     77    !$OMP END PARALLEL
     78
     79  END SUBROUTINE iniphysiq
    8380
    8481END MODULE iniphysiq_mod

  • LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/phylmd/ce0l.F90

    r5117 r5118  
    11PROGRAM ce0l
    22
    3 !-------------------------------------------------------------------------------
    4 ! Purpose: Initial states and boundary conditions files creation:
    5 !     * start.nc    for dynamics    (using etat0dyn     routine)
    6 !     * startphy.nc for physics     (using etat0phys    routine)
    7 !     * limit.nc    for forced runs (using limit_netcdf routine)
    8 !-------------------------------------------------------------------------------
    9 ! Notes:
    10 !     * extrap=.T. (default) for data extrapolation, like for the SSTs when file
    11 !                   does contain ocean points only.
    12 !     * "masque" can be:
    13 !       - read from file "o2a.nc"          (for coupled runs).
    14 !       - read from file "startphy0.nc"    (from a previous run).
    15 !       - created in etat0phys or etat0dyn (for forced  runs).
    16 !     It is then passed to limit_netcdf to ensure consistancy.
    17 !-------------------------------------------------------------------------------
     3  !-------------------------------------------------------------------------------
     4  ! Purpose: Initial states and boundary conditions files creation:
     5  !     * start.nc    for dynamics    (using etat0dyn     routine)
     6  !     * startphy.nc for physics     (using etat0phys    routine)
     7  !     * limit.nc    for forced runs (using limit_netcdf routine)
     8  !-------------------------------------------------------------------------------
     9  ! Notes:
     10  !     * extrap=.T. (default) for data extrapolation, like for the SSTs when file
     11  !                   does contain ocean points only.
     12  !     * "masque" can be:
     13  !       - read from file "o2a.nc"          (for coupled runs).
     14  !       - read from file "startphy0.nc"    (from a previous run).
     15  !       - created in etat0phys or etat0dyn (for forced  runs).
     16  !     It is then passed to limit_netcdf to ensure consistancy.
     17  !-------------------------------------------------------------------------------
    1818  USE ioipsl, ONLY: ioconf_calendar, getin, flininfo, flinopen, flinget, flinclo
    19   USE control_mod,    ONLY: day_step, dayref, nsplit_phys
    20   USE etat0dyn,       ONLY: etat0dyn_netcdf
    21   USE etat0phys,      ONLY: etat0phys_netcdf
    22   USE limit,          ONLY: limit_netcdf
    23   USE netcdf,         ONLY: nf90_open, nf90_nowrite, nf90_close, nf90_noerr,    &
    24          nf90_inquire_dimension, nf90_inq_dimid, nf90_inq_varid, nf90_get_var
    25   USE infotrac,       ONLY: init_infotrac
    26   USE dimphy,         ONLY: klon
     19  USE control_mod, ONLY: day_step, dayref, nsplit_phys
     20  USE etat0dyn, ONLY: etat0dyn_netcdf
     21  USE etat0phys, ONLY: etat0phys_netcdf
     22  USE limit, ONLY: limit_netcdf
     23  USE netcdf, ONLY: nf90_open, nf90_nowrite, nf90_close, nf90_noerr, &
     24          nf90_inquire_dimension, nf90_inq_dimid, nf90_inq_varid, nf90_get_var
     25  USE infotrac, ONLY: init_infotrac
     26  USE dimphy, ONLY: klon
    2727  USE test_disvert_m, ONLY: test_disvert
    28   USE lmdz_filtreg,    ONLY: inifilr
    29   USE iniphysiq_mod,  ONLY: iniphysiq
    30   USE mod_const_mpi,  ONLY: comm_lmdz
     28  USE lmdz_filtreg, ONLY: inifilr
     29  USE iniphysiq_mod, ONLY: iniphysiq
     30  USE mod_const_mpi, ONLY: comm_lmdz
    3131
    3232#ifdef CPP_PARA
    … …  
    4040
    4141  USE comconst_mod, ONLY: cpp, daysec, dtphys, dtvr, g, kappa, omeg, r, rad, &
    42                           pi, jmp1
     42          pi, jmp1
    4343  USE logic_mod, ONLY: iflag_phys, ok_etat0, ok_limit
    4444  USE comvert_mod, ONLY: pa, preff, pressure_exner
    4545  USE temps_mod, ONLY: calend, day_ini, dt
    4646  USE lmdz_mpi
     47  USE lmdz_iniprint, ONLY: lunout, prt_level
    4748
    4849  IMPLICIT NONE
    4950
    50 !-------------------------------------------------------------------------------
    51 ! Local variables:
     51  !-------------------------------------------------------------------------------
     52  ! Local variables:
    5253  include "dimensions.h"
    5354  include "paramet.h"
    5455  include "comgeom2.h"
    55   include "iniprint.h"
    56  
    57   REAL               :: masque(iip1,jjp1)             !--- CONTINENTAL MASK
    58   REAL               :: phis  (iip1,jjp1)             !--- GROUND GEOPOTENTIAL
    59   CHARACTER(LEN=256) :: modname, fmt, calnd           !--- CALENDAR TYPE
    60   LOGICAL            :: use_filtre_fft
    61   LOGICAL, PARAMETER :: extrap=.FALSE.
    62 
    63 !--- Local variables for ocean mask reading:
    64   INTEGER            :: nid_o2a, iml_omask, jml_omask, j
    65   INTEGER            :: fid, iret, llm_tmp, ttm_tmp, itaul(1)
    66   REAL, ALLOCATABLE  :: lon_omask(:,:), dlon_omask(:), ocemask(:,:)
    67   REAL, ALLOCATABLE  :: lat_omask(:,:), dlat_omask(:), ocetmp (:,:)
    68   REAL               :: date, lev(1)
    69 
    70 !--- Local variables for land mask from startphy0 file reading
    71   INTEGER            :: nid_sta, nid_nph, nid_msk, nphys
    72   REAL, ALLOCATABLE  :: masktmp(:)
     56
     57  REAL :: masque(iip1, jjp1)             !--- CONTINENTAL MASK
     58  REAL :: phis  (iip1, jjp1)             !--- GROUND GEOPOTENTIAL
     59  CHARACTER(LEN = 256) :: modname, fmt, calnd           !--- CALENDAR TYPE
     60  LOGICAL :: use_filtre_fft
     61  LOGICAL, PARAMETER :: extrap = .FALSE.
     62
     63  !--- Local variables for ocean mask reading:
     64  INTEGER :: nid_o2a, iml_omask, jml_omask, j
     65  INTEGER :: fid, iret, llm_tmp, ttm_tmp, itaul(1)
     66  REAL, ALLOCATABLE :: lon_omask(:, :), dlon_omask(:), ocemask(:, :)
     67  REAL, ALLOCATABLE :: lat_omask(:, :), dlat_omask(:), ocetmp (:, :)
     68  REAL :: date, lev(1)
     69
     70  !--- Local variables for land mask from startphy0 file reading
     71  INTEGER :: nid_sta, nid_nph, nid_msk, nphys
     72  REAL, ALLOCATABLE :: masktmp(:)
    7373
    7474#ifdef CPP_PARA
    7575  INTEGER ierr
    7676#else
    77 ! for iniphysiq in serial mode
    78   INTEGER,PARAMETER :: mpi_rank=0
    79   INTEGER :: distrib_phys(mpi_rank:mpi_rank)=(jjm-1)*iim+2
    80 #endif
    81 !-------------------------------------------------------------------------------
    82   modname="ce0l"
    83 
    84 !--- Constants
    85   pi     = 4. * ATAN(1.)
    86   rad    = 6371229.
     77  ! for iniphysiq in serial mode
     78  INTEGER, PARAMETER :: mpi_rank = 0
     79  INTEGER :: distrib_phys(mpi_rank:mpi_rank) = (jjm - 1) * iim + 2
     80#endif
     81  !-------------------------------------------------------------------------------
     82  modname = "ce0l"
     83
     84  !--- Constants
     85  pi = 4. * ATAN(1.)
     86  rad = 6371229.
    8787  daysec = 86400.
    88   omeg   = 2.*pi/daysec
    89   g      = 9.8
    90   kappa  = 0.2857143
    91   cpp    = 1004.70885
    92   jmp1   = jjm + 1
    93   preff   = 101325.
    94   pa      = 50000.
    95 
    96   CALL conf_gcm( 99, .TRUE. )
    97   dtvr = daysec/REAL(day_step)
    98   WRITE(lunout,*)'dtvr',dtvr
     88  omeg = 2. * pi / daysec
     89  g = 9.8
     90  kappa = 0.2857143
     91  cpp = 1004.70885
     92  jmp1 = jjm + 1
     93  preff = 101325.
     94  pa = 50000.
     95
     96  CALL conf_gcm(99, .TRUE.)
     97  dtvr = daysec / REAL(day_step)
     98  WRITE(lunout, *)'dtvr', dtvr
    9999  CALL iniconst()
    100100  CALL inigeom()
    101101
    102 !--- Calendar choice
    103   calnd='gregorian'
     102  !--- Calendar choice
     103  calnd = 'gregorian'
    104104  SELECT CASE(calend)
    105     CASE('earth_360d');CALL ioconf_calendar('360_day');   calnd='with 360 days/year'
    106     CASE('earth_365d');CALL ioconf_calendar('noleap'); calnd='with no leap year'
    107     CASE('earth_366d');CALL ioconf_calendar('366d');   calnd='with leap years only'
    108     CASE('gregorian'); CALL ioconf_calendar('gregorian')
    109     CASE('standard');  CALL ioconf_calendar('gregorian')
    110     CASE('julian');    CALL ioconf_calendar('julian'); calnd='julian'
    111     CASE('proleptic_gregorian'); CALL ioconf_calendar('gregorian')
     105  CASE('earth_360d');CALL ioconf_calendar('360_day');   calnd = 'with 360 days/year'
     106  CASE('earth_365d');CALL ioconf_calendar('noleap'); calnd = 'with no leap year'
     107  CASE('earth_366d');CALL ioconf_calendar('366d');   calnd = 'with leap years only'
     108  CASE('gregorian'); CALL ioconf_calendar('gregorian')
     109  CASE('standard');  CALL ioconf_calendar('gregorian')
     110  CASE('julian');    CALL ioconf_calendar('julian'); calnd = 'julian'
     111  CASE('proleptic_gregorian'); CALL ioconf_calendar('gregorian')
    112112  !--- DC Bof...  => IOIPSL a mettre a jour: proleptic_gregorian /= gregorian
    113     CASE DEFAULT
    114       CALL abort_gcm('ce0l','Bad choice for calendar',1)
     113  CASE DEFAULT
     114    CALL abort_gcm('ce0l', 'Bad choice for calendar', 1)
    115115  END SELECT
    116   WRITE(lunout,*)'CHOSEN CALENDAR: Earth '//TRIM(calnd)
     116  WRITE(lunout, *)'CHOSEN CALENDAR: Earth ' // TRIM(calnd)
    117117
    118118#ifdef CPP_PARA
    … …  
    123123  CALL init_mod_hallo()
    124124#endif
    125   WRITE(lunout,*)'---> klon=',klon
    126 
    127 !--- Tracers initializations
     125  WRITE(lunout, *)'---> klon=', klon
     126
     127  !--- Tracers initializations
    128128  CALL init_infotrac()
    129129
    130130  CALL inifilr()
    131   CALL iniphysiq(iim,jjm,llm, &
    132                  distrib_phys(mpi_rank),comm_lmdz, &
    133                  daysec,day_ini,dtphys/nsplit_phys, &
    134                  rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys)
     131  CALL iniphysiq(iim, jjm, llm, &
     132          distrib_phys(mpi_rank), comm_lmdz, &
     133          daysec, day_ini, dtphys / nsplit_phys, &
     134          rlatu, rlatv, rlonu, rlonv, aire, cu, cv, rad, g, r, cpp, iflag_phys)
    135135  IF(pressure_exner) CALL test_disvert
    136136
    … …  
    138138  IF (mpi_rank==0.AND.omp_rank==0) THEN
    139139#endif
    140   use_filtre_fft=.FALSE.
    141   CALL getin('use_filtre_fft',use_filtre_fft)
     140  use_filtre_fft = .FALSE.
     141  CALL getin('use_filtre_fft', use_filtre_fft)
    142142  IF(use_filtre_fft) THEN
    143      WRITE(lunout,*)"FFT filter not available for sequential dynamics."
    144      WRITE(lunout,*)"Your setting of variable use_filtre_fft is not used."
     143    WRITE(lunout, *)"FFT filter not available for sequential dynamics."
     144    WRITE(lunout, *)"Your setting of variable use_filtre_fft is not used."
    145145  ENDIF
    146146
    147 !--- LAND MASK. THREE CASES:
    148 !   1) read from ocean model    file "o2a.nc"    (coupled runs)
    149 !   2) read from previous run   file="startphy0.nc"
    150 !   3) computed from topography file "Relief.nc" (masque(:,:)=-99999.)
    151 ! In the first case, the mask from the ocean model is used compute the
    152 ! weights to ensure ocean fractions are the same for atmosphere and ocean.
    153 !*******************************************************************************
     147  !--- LAND MASK. THREE CASES:
     148  !   1) read from ocean model    file "o2a.nc"    (coupled runs)
     149  !   2) read from previous run   file="startphy0.nc"
     150  !   3) computed from topography file "Relief.nc" (masque(:,:)=-99999.)
     151  ! In the first case, the mask from the ocean model is used compute the
     152  ! weights to ensure ocean fractions are the same for atmosphere and ocean.
     153  !*******************************************************************************
    154154  IF(nf90_open("o2a.nc", nf90_nowrite, nid_o2a)==nf90_noerr) THEN
    155     iret=nf90_close(nid_o2a)
    156     WRITE(lunout,*)'BEWARE !! Ocean mask "o2a.nc" file found'
    157     WRITE(lunout,*)'Coupled run.'
     155    iret = nf90_close(nid_o2a)
     156    WRITE(lunout, *)'BEWARE !! Ocean mask "o2a.nc" file found'
     157    WRITE(lunout, *)'Coupled run.'
    158158    CALL flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp, nid_o2a)
    159159    IF(iml_omask/=iim .OR.jml_omask/=jjp1) THEN
    160       WRITE(lunout,*)'Mismatching dimensions for ocean mask'
    161       WRITE(lunout,*)'iim  = ',iim ,' iml_omask = ',iml_omask
    162       WRITE(lunout,*)'jjp1 = ',jjp1,' jml_omask = ',jml_omask
    163       CALL abort_gcm(modname,'',1)
    164     END IF
    165     ALLOCATE(ocemask(iim,jjp1),lon_omask(iim,jjp1),dlon_omask(iim ))
    166     ALLOCATE(ocetmp (iim,jjp1),lat_omask(iim,jjp1),dlat_omask(jjp1))
    167     CALL flinopen("o2a.nc", .FALSE.,iml_omask,jml_omask,llm_tmp,              &
    168                   lon_omask,lat_omask,lev,ttm_tmp,itaul,date,dt,fid)
    169     CALL flinget(fid, "OceMask",    iim,jjp1,llm_tmp,ttm_tmp,1,1,ocetmp)
     160      WRITE(lunout, *)'Mismatching dimensions for ocean mask'
     161      WRITE(lunout, *)'iim  = ', iim, ' iml_omask = ', iml_omask
     162      WRITE(lunout, *)'jjp1 = ', jjp1, ' jml_omask = ', jml_omask
     163      CALL abort_gcm(modname, '', 1)
     164    END IF
     165    ALLOCATE(ocemask(iim, jjp1), lon_omask(iim, jjp1), dlon_omask(iim))
     166    ALLOCATE(ocetmp (iim, jjp1), lat_omask(iim, jjp1), dlat_omask(jjp1))
     167    CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp, &
     168            lon_omask, lat_omask, lev, ttm_tmp, itaul, date, dt, fid)
     169    CALL flinget(fid, "OceMask", iim, jjp1, llm_tmp, ttm_tmp, 1, 1, ocetmp)
    170170    CALL flinclo(fid)
    171     dlon_omask(1:iim ) = lon_omask(1:iim,1)
    172     dlat_omask(1:jjp1) = lat_omask(1,1:jjp1)
     171    dlon_omask(1:iim) = lon_omask(1:iim, 1)
     172    dlat_omask(1:jjp1) = lat_omask(1, 1:jjp1)
    173173    ocemask = ocetmp
    174174    IF(dlat_omask(1)<dlat_omask(jml_omask)) THEN
    175        DO j=1,jjp1
    176           ocemask(:,j) = ocetmp(:,jjp1-j+1)
    177        END DO
    178     END IF
    179     DEALLOCATE(ocetmp,lon_omask,lat_omask,dlon_omask,dlat_omask)
     175      DO j = 1, jjp1
     176        ocemask(:, j) = ocetmp(:, jjp1 - j + 1)
     177      END DO
     178    END IF
     179    DEALLOCATE(ocetmp, lon_omask, lat_omask, dlon_omask, dlat_omask)
    180180    IF(prt_level>=1) THEN
    181       WRITE(fmt,"(i4,'i1)')")iim ; fmt='('//ADJUSTL(fmt)
    182       WRITE(lunout,*)'OCEAN MASK :'
    183       WRITE(lunout,fmt) NINT(ocemask)
    184     END IF
    185     masque(1:iim,:)=1.-ocemask(:,:)
    186     masque(iip1 ,:)=masque(1,:)
     181      WRITE(fmt, "(i4,'i1)')")iim ; fmt = '(' // ADJUSTL(fmt)
     182      WRITE(lunout, *)'OCEAN MASK :'
     183      WRITE(lunout, fmt) NINT(ocemask)
     184    END IF
     185    masque(1:iim, :) = 1. - ocemask(:, :)
     186    masque(iip1, :) = masque(1, :)
    187187    DEALLOCATE(ocemask)
    188188  ELSE IF(nf90_open("startphy0.nc", nf90_nowrite, nid_sta)==nf90_noerr) THEN
    189     WRITE(lunout,*)'BEWARE !! File "startphy0.nc" found.'
    190     WRITE(lunout,*)'Getting the land mask from a previous run.'
    191     iret=nf90_inq_dimid(nid_sta,'points_physiques',nid_nph)
    192     iret=nf90_inquire_dimension(nid_sta,nid_nph,len=nphys)
     189    WRITE(lunout, *)'BEWARE !! File "startphy0.nc" found.'
     190    WRITE(lunout, *)'Getting the land mask from a previous run.'
     191    iret = nf90_inq_dimid(nid_sta, 'points_physiques', nid_nph)
     192    iret = nf90_inquire_dimension(nid_sta, nid_nph, len = nphys)
    193193    IF(nphys/=klon) THEN
    194       WRITE(lunout,*)'Mismatching dimensions for land mask'
    195       WRITE(lunout,*)'nphys  = ',nphys ,' klon = ',klon
    196       iret=nf90_close(nid_sta)
    197       CALL abort_gcm(modname,'',1)
     194      WRITE(lunout, *)'Mismatching dimensions for land mask'
     195      WRITE(lunout, *)'nphys  = ', nphys, ' klon = ', klon
     196      iret = nf90_close(nid_sta)
     197      CALL abort_gcm(modname, '', 1)
    198198    END IF
    199199    ALLOCATE(masktmp(klon))
    200     iret=nf90_inq_varid(nid_sta,'masque',nid_msk)
    201     iret=nf90_get_var(nid_sta,nid_msk,masktmp)
    202     iret=nf90_close(nid_sta)
    203     CALL gr_fi_dyn(1,klon,iip1,jjp1,masktmp,masque)
     200    iret = nf90_inq_varid(nid_sta, 'masque', nid_msk)
     201    iret = nf90_get_var(nid_sta, nid_msk, masktmp)
     202    iret = nf90_close(nid_sta)
     203    CALL gr_fi_dyn(1, klon, iip1, jjp1, masktmp, masque)
    204204    IF(prt_level>=1) THEN
    205       WRITE(fmt,"(i4,'i1)')")iip1 ; fmt='('//ADJUSTL(fmt)
    206       WRITE(lunout,*)'LAND MASK :'
    207       WRITE(lunout,fmt) NINT(masque)
     205      WRITE(fmt, "(i4,'i1)')")iip1 ; fmt = '(' // ADJUSTL(fmt)
     206      WRITE(lunout, *)'LAND MASK :'
     207      WRITE(lunout, fmt) NINT(masque)
    208208    END IF
    209209    DEALLOCATE(masktmp)
    210210  ELSE
    211     WRITE(lunout,*)'BEWARE !! No ocean mask "o2a.nc" file or "startphy0.nc" file found'
    212     WRITE(lunout,*)'Land mask will be built from the topography file.'
    213     masque(:,:)=-99999.
    214   END IF
    215   phis(:,:)=-99999.
     211    WRITE(lunout, *)'BEWARE !! No ocean mask "o2a.nc" file or "startphy0.nc" file found'
     212    WRITE(lunout, *)'Land mask will be built from the topography file.'
     213    masque(:, :) = -99999.
     214  END IF
     215  phis(:, :) = -99999.
    216216
    217217  IF(ok_etat0) THEN
    218     WRITE(lunout,'(//)')
    219     WRITE(lunout,*) '  ************************  '
    220     WRITE(lunout,*) '  ***  etat0phy_netcdf ***  '
    221     WRITE(lunout,*) '  ************************  '
    222     CALL etat0phys_netcdf(masque,phis)
    223     WRITE(lunout,'(//)')
    224     WRITE(lunout,*) '  ************************  '
    225     WRITE(lunout,*) '  ***  etat0dyn_netcdf ***  '
    226     WRITE(lunout,*) '  ************************  '
    227     CALL etat0dyn_netcdf(masque,phis)
     218    WRITE(lunout, '(//)')
     219    WRITE(lunout, *) '  ************************  '
     220    WRITE(lunout, *) '  ***  etat0phy_netcdf ***  '
     221    WRITE(lunout, *) '  ************************  '
     222    CALL etat0phys_netcdf(masque, phis)
     223    WRITE(lunout, '(//)')
     224    WRITE(lunout, *) '  ************************  '
     225    WRITE(lunout, *) '  ***  etat0dyn_netcdf ***  '
     226    WRITE(lunout, *) '  ************************  '
     227    CALL etat0dyn_netcdf(masque, phis)
    228228  END IF
    229229
    230230  IF(ok_limit) THEN
    231     WRITE(lunout,'(//)')
    232     WRITE(lunout,*) '  *********************  '
    233     WRITE(lunout,*) '  ***  Limit_netcdf ***  '
    234     WRITE(lunout,*) '  *********************  '
    235     WRITE(lunout,'(//)')
    236     CALL limit_netcdf(masque,phis,extrap)
    237   END IF
    238 
    239   WRITE(lunout,'(//)')
    240   WRITE(lunout,*) '  ***************************  '
    241   WRITE(lunout,*) '  ***  grilles_gcm_netcdf ***  '
    242   WRITE(lunout,*) '  ***************************  '
    243   WRITE(lunout,'(//)')
    244   CALL grilles_gcm_netcdf_sub(masque,phis)
     231    WRITE(lunout, '(//)')
     232    WRITE(lunout, *) '  *********************  '
     233    WRITE(lunout, *) '  ***  Limit_netcdf ***  '
     234    WRITE(lunout, *) '  *********************  '
     235    WRITE(lunout, '(//)')
     236    CALL limit_netcdf(masque, phis, extrap)
     237  END IF
     238
     239  WRITE(lunout, '(//)')
     240  WRITE(lunout, *) '  ***************************  '
     241  WRITE(lunout, *) '  ***  grilles_gcm_netcdf ***  '
     242  WRITE(lunout, *) '  ***************************  '
     243  WRITE(lunout, '(//)')
     244  CALL grilles_gcm_netcdf_sub(masque, phis)
    245245
    246246#ifdef CPP_PARA

  • LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/phylmd/etat0dyn_netcdf.F90

    r5117 r5118  
    4040  USE temps_mod, ONLY: annee_ref, day_ref, itau_dyn, itau_phy, start_time
    4141  USE lmdz_strings, ONLY: strLower
     42  USE lmdz_iniprint, ONLY: lunout, prt_level
    4243
    4344  IMPLICIT NONE
    … …  
    4647  PUBLIC :: etat0dyn_netcdf
    4748
    48   include "iniprint.h"
    4949  include "dimensions.h"
    5050  include "paramet.h"

  • LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90

    r5117 r5118  
    1 
    21! $Id$
    32
    43MODULE etat0phys
    54
    6 !*******************************************************************************
    7 ! Purpose: Create physical initial state using atmospheric fields from a
    8 !          database of atmospheric to initialize the model.
    9 !-------------------------------------------------------------------------------
    10 ! Comments:
    11 
    12 !    *  This module is designed to work for Earth (and with ioipsl)
    13 
    14 !    *  etat0phys_netcdf routine can access to NetCDF data through subroutines:
    15 !         "start_init_phys" for variables contained in file "ECPHY.nc":
    16 !            'ST'     : Surface temperature
    17 !            'CDSW'   : Soil moisture
    18 !         "start_init_orog" for variables contained in file "Relief.nc":
    19 !            'RELIEF' : High resolution orography
    20 
    21 !    * The land mask and corresponding weights can be:
    22 !      1) computed using the ocean mask from the ocean model (to ensure ocean
    23 !         fractions are the same for atmosphere and ocean) for coupled runs.
    24 !         File name: "o2a.nc"  ;  variable name: "OceMask"
    25 !      2) computed from topography file "Relief.nc" for forced runs.
    26 
    27 !    * Allowed values for read_climoz flag are 0, 1 and 2:
    28 !      0: do not read an ozone climatology
    29 !      1: read a single ozone climatology that will be used day and night
    30 !      2: read two ozone climatologies, the average day and night climatology
    31 !         and the daylight climatology
    32 !-------------------------------------------------------------------------------
    33 !    * There is a big mess with the longitude size. Should it be iml or iml+1 ?
    34 !  I have chosen to use the iml+1 as an argument to this routine and we declare
    35 !  internaly smaller fields when needed. This needs to be cleared once and for
    36 !  all in LMDZ. A convention is required.
    37 !-------------------------------------------------------------------------------
    38 
    39   USE ioipsl,             ONLY: flininfo, flinopen, flinget, flinclo
    40   USE lmdz_assert_eq,        ONLY: assert_eq
    41   USE dimphy,             ONLY: klon
    42   USE conf_dat_m,         ONLY: conf_dat2d
     5  !*******************************************************************************
     6  ! Purpose: Create physical initial state using atmospheric fields from a
     7  !          database of atmospheric to initialize the model.
     8  !-------------------------------------------------------------------------------
     9  ! Comments:
     10
     11  !    *  This module is designed to work for Earth (and with ioipsl)
     12
     13  !    *  etat0phys_netcdf routine can access to NetCDF data through subroutines:
     14  !         "start_init_phys" for variables contained in file "ECPHY.nc":
     15  !            'ST'     : Surface temperature
     16  !            'CDSW'   : Soil moisture
     17  !         "start_init_orog" for variables contained in file "Relief.nc":
     18  !            'RELIEF' : High resolution orography
     19
     20  !    * The land mask and corresponding weights can be:
     21  !      1) computed using the ocean mask from the ocean model (to ensure ocean
     22  !         fractions are the same for atmosphere and ocean) for coupled runs.
     23  !         File name: "o2a.nc"  ;  variable name: "OceMask"
     24  !      2) computed from topography file "Relief.nc" for forced runs.
     25
     26  !    * Allowed values for read_climoz flag are 0, 1 and 2:
     27  !      0: do not read an ozone climatology
     28  !      1: read a single ozone climatology that will be used day and night
     29  !      2: read two ozone climatologies, the average day and night climatology
     30  !         and the daylight climatology
     31  !-------------------------------------------------------------------------------
     32  !    * There is a big mess with the longitude size. Should it be iml or iml+1 ?
     33  !  I have chosen to use the iml+1 as an argument to this routine and we declare
     34  !  internaly smaller fields when needed. This needs to be cleared once and for
     35  !  all in LMDZ. A convention is required.
     36  !-------------------------------------------------------------------------------
     37
     38  USE ioipsl, ONLY: flininfo, flinopen, flinget, flinclo
     39  USE lmdz_assert_eq, ONLY: assert_eq
     40  USE dimphy, ONLY: klon
     41  USE conf_dat_m, ONLY: conf_dat2d
    4342  USE phys_state_var_mod, ONLY: zmea, zstd, zsig, zgam, zthe, zpic, zval, z0m, &
    44           solsw, solswfdiff, radsol, t_ancien, wake_deltat, wake_s,  rain_fall, qsol, z0h, &
    45           sollw,sollwdown, rugoro, q_ancien, wake_deltaq, wake_pe, snow_fall, ratqs,w01, &
    46     sig1, ftsol, clwcon, fm_therm, wake_Cstar,  pctsrf,  entr_therm,radpas, f0,&
    47     zmax0,fevap, rnebcon,falb_dir, falb_dif, wake_fip,    agesno,  detr_therm, pbl_tke, &
    48     phys_state_var_init, ql_ancien, qs_ancien, prlw_ancien, prsw_ancien, &
    49     prw_ancien, u10m,v10m, treedrg, u_ancien, v_ancien, wake_delta_pbl_TKE, wake_dens, &
    50     ale_bl, ale_bl_trig, alp_bl, &
    51     ale_wake, ale_bl_stat, AWAKE_S
     43          solsw, solswfdiff, radsol, t_ancien, wake_deltat, wake_s, rain_fall, qsol, z0h, &
     44          sollw, sollwdown, rugoro, q_ancien, wake_deltaq, wake_pe, snow_fall, ratqs, w01, &
     45          sig1, ftsol, clwcon, fm_therm, wake_Cstar, pctsrf, entr_therm, radpas, f0, &
     46          zmax0, fevap, rnebcon, falb_dir, falb_dif, wake_fip, agesno, detr_therm, pbl_tke, &
     47          phys_state_var_init, ql_ancien, qs_ancien, prlw_ancien, prsw_ancien, &
     48          prw_ancien, u10m, v10m, treedrg, u_ancien, v_ancien, wake_delta_pbl_TKE, wake_dens, &
     49          ale_bl, ale_bl_trig, alp_bl, &
     50          ale_wake, ale_bl_stat, AWAKE_S
    5251
    5352  USE comconst_mod, ONLY: pi, dtvr
     53  USE lmdz_iniprint, ONLY: lunout, prt_level
    5454
    5555  PRIVATE
    5656  PUBLIC :: etat0phys_netcdf
    5757
    58   include "iniprint.h"
    5958  include "dimensions.h"
    6059  include "paramet.h"
    … …  
    6463  REAL, SAVE :: deg2rad
    6564  REAL, SAVE, ALLOCATABLE :: tsol(:)
    66   INTEGER,            SAVE      :: iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys
    67   REAL, ALLOCATABLE,  SAVE      :: lon_phys(:,:), lat_phys(:,:), levphys_ini(:)
    68   CHARACTER(LEN=256), PARAMETER :: oroparam="oro_params.nc"
    69   CHARACTER(LEN=256), PARAMETER :: orofname="Relief.nc", orogvar="RELIEF"
    70   CHARACTER(LEN=256), PARAMETER :: phyfname="ECPHY.nc",  psrfvar="SP"
    71   CHARACTER(LEN=256), PARAMETER :: qsolvar="CDSW",       tsrfvar="ST"
     65  INTEGER, SAVE :: iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys
     66  REAL, ALLOCATABLE, SAVE :: lon_phys(:, :), lat_phys(:, :), levphys_ini(:)
     67  CHARACTER(LEN = 256), PARAMETER :: oroparam = "oro_params.nc"
     68  CHARACTER(LEN = 256), PARAMETER :: orofname = "Relief.nc", orogvar = "RELIEF"
     69  CHARACTER(LEN = 256), PARAMETER :: phyfname = "ECPHY.nc", psrfvar = "SP"
     70  CHARACTER(LEN = 256), PARAMETER :: qsolvar = "CDSW", tsrfvar = "ST"
    7271
    7372
    … …  
    7574
    7675
    77 !-------------------------------------------------------------------------------
    78 
    79 SUBROUTINE etat0phys_netcdf(masque, phis)
    80 
    81 !-------------------------------------------------------------------------------
    82 ! Purpose: Creates initial states
    83 !-------------------------------------------------------------------------------
    84 ! Notes:  1) This routine is designed to work for Earth
    85 !         2) If masque(:,:)/=-99999., masque and phis are already known.
    86 !         Otherwise: compute it.
    87 !-------------------------------------------------------------------------------
    88   USE control_mod
    89   USE fonte_neige_mod
    90   USE pbl_surface_mod
    91   USE regr_horiz_time_climoz_m, ONLY: regr_horiz_time_climoz
    92   USE indice_sol_mod
    93   USE conf_phys_m, ONLY: conf_phys
    94   USE init_ssrf_m, ONLY: start_init_subsurf
    95   USE phys_state_var_mod, ONLY: beta_aridity, delta_tsurf, awake_dens, cv_gen, &
    96        ratqs_inter_, rneb_ancien
    97   !use ioipsl_getincom
    98   IMPLICIT NONE
    99 !-------------------------------------------------------------------------------
    100 ! Arguments:
    101   REAL,    INTENT(INOUT) :: masque(:,:) !--- Land mask           dim(iip1,jjp1)
    102   REAL,    INTENT(INOUT) :: phis  (:,:) !--- Ground geopotential dim(iip1,jjp1)
    103 !-------------------------------------------------------------------------------
    104 ! Local variables:
    105   CHARACTER(LEN=256) :: modname="etat0phys_netcdf", fmt
    106   INTEGER            :: i, j, l, ji, iml, jml
    107   LOGICAL            :: read_mask
    108   REAL               :: phystep, dummy
    109   REAL, DIMENSION(SIZE(masque,1),SIZE(masque,2)) :: masque_tmp,phiso
    110   REAL, DIMENSION(klon)               :: sn, rugmer, run_off_lic_0, fder
    111   REAL, DIMENSION(klon,nbsrf)         :: qsurf, snsrf
    112   REAL, DIMENSION(klon,nsoilmx,nbsrf) :: tsoil
    113 
    114 !--- Arguments for conf_phys
    115   LOGICAL :: ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, callstats
    116   REAL    :: solarlong0, seuil_inversion, fact_cldcon, facttemps
    117   LOGICAL :: ok_newmicro
    118   INTEGER :: iflag_radia, iflag_cldcon, iflag_ratqs
    119   REAL    :: ratqsbas, ratqshaut, tau_ratqs
    120   LOGICAL :: ok_ade, ok_aie, ok_volcan, ok_alw, ok_cdnc, aerosol_couple, chemistry_couple
    121   INTEGER :: flag_aerosol
    122   INTEGER :: flag_aerosol_strat
    123   INTEGER :: flag_volc_surfstrat
    124   LOGICAL :: flag_aer_feedback
    125   LOGICAL :: flag_bc_internal_mixture
    126   REAL    :: bl95_b0, bl95_b1
    127   INTEGER :: read_climoz                        !--- Read ozone climatology
    128   REAL    :: alp_offset
    129   LOGICAL :: filtre_oro=.FALSE.
    130 
    131   INCLUDE "compbl.h"
    132   INCLUDE "alpale.h"
    133  
    134   deg2rad= pi/180.0
    135   iml=assert_eq(SIZE(masque,1),SIZE(phis,1),TRIM(modname)//" iml")
    136   jml=assert_eq(SIZE(masque,2),SIZE(phis,2),TRIM(modname)//" jml")
    137 
    138 ! Physics configuration
    139 !*******************************************************************************
    140   CALL conf_phys(  ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES,     &
    141                    callstats,                                           &
    142                    solarlong0,seuil_inversion,                          &
    143                    fact_cldcon, facttemps,ok_newmicro,iflag_radia,      &
    144                    iflag_cldcon,                                        &
    145                    ratqsbas,ratqshaut,tau_ratqs,            &
    146                    ok_ade, ok_aie, ok_alw, ok_cdnc, ok_volcan, flag_volc_surfstrat,     &
    147                    aerosol_couple, chemistry_couple, flag_aerosol, flag_aerosol_strat,  &
    148                    flag_aer_feedback, flag_bc_internal_mixture, bl95_b0, bl95_b1,       &
    149                    read_climoz, alp_offset)
    150   CALL phys_state_var_init(read_climoz)
    151 
    152 !--- Initial atmospheric CO2 conc. from .def file
    153   co2_ppm0 = co2_ppm
    154 
    155 ! Compute ground geopotential, sub-cells quantities and possibly the mask.
    156 !*******************************************************************************
    157   read_mask=ANY(masque/=-99999.); masque_tmp=masque
    158   CALL start_init_orog(rlonv, rlatu, phis, masque_tmp)
    159 
    160   CALL getin('filtre_oro',filtre_oro)
    161   IF (filtre_oro) CALL filtreoro(size(phis,1),size(phis,2),phis,masque_tmp,rlatu)
    162 
    163   WRITE(fmt,"(i4,'i1)')")iml ; fmt='('//ADJUSTL(fmt)
    164   IF(.NOT.read_mask) THEN                       !--- Keep mask form orography
    165     masque=masque_tmp
    166     IF(prt_level>=1) THEN
    167       WRITE(lunout,*)'BUILT MASK :'
    168       WRITE(lunout,fmt) NINT(masque)
     76  !-------------------------------------------------------------------------------
     77
     78  SUBROUTINE etat0phys_netcdf(masque, phis)
     79
     80    !-------------------------------------------------------------------------------
     81    ! Purpose: Creates initial states
     82    !-------------------------------------------------------------------------------
     83    ! Notes:  1) This routine is designed to work for Earth
     84    !         2) If masque(:,:)/=-99999., masque and phis are already known.
     85    !         Otherwise: compute it.
     86    !-------------------------------------------------------------------------------
     87    USE control_mod
     88    USE fonte_neige_mod
     89    USE pbl_surface_mod
     90    USE regr_horiz_time_climoz_m, ONLY: regr_horiz_time_climoz
     91    USE indice_sol_mod
     92    USE conf_phys_m, ONLY: conf_phys
     93    USE init_ssrf_m, ONLY: start_init_subsurf
     94    USE phys_state_var_mod, ONLY: beta_aridity, delta_tsurf, awake_dens, cv_gen, &
     95            ratqs_inter_, rneb_ancien
     96    !use ioipsl_getincom
     97    IMPLICIT NONE
     98    !-------------------------------------------------------------------------------
     99    ! Arguments:
     100    REAL, INTENT(INOUT) :: masque(:, :) !--- Land mask           dim(iip1,jjp1)
     101    REAL, INTENT(INOUT) :: phis  (:, :) !--- Ground geopotential dim(iip1,jjp1)
     102    !-------------------------------------------------------------------------------
     103    ! Local variables:
     104    CHARACTER(LEN = 256) :: modname = "etat0phys_netcdf", fmt
     105    INTEGER :: i, j, l, ji, iml, jml
     106    LOGICAL :: read_mask
     107    REAL :: phystep, dummy
     108    REAL, DIMENSION(SIZE(masque, 1), SIZE(masque, 2)) :: masque_tmp, phiso
     109    REAL, DIMENSION(klon) :: sn, rugmer, run_off_lic_0, fder
     110    REAL, DIMENSION(klon, nbsrf) :: qsurf, snsrf
     111    REAL, DIMENSION(klon, nsoilmx, nbsrf) :: tsoil
     112
     113    !--- Arguments for conf_phys
     114    LOGICAL :: ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, callstats
     115    REAL :: solarlong0, seuil_inversion, fact_cldcon, facttemps
     116    LOGICAL :: ok_newmicro
     117    INTEGER :: iflag_radia, iflag_cldcon, iflag_ratqs
     118    REAL :: ratqsbas, ratqshaut, tau_ratqs
     119    LOGICAL :: ok_ade, ok_aie, ok_volcan, ok_alw, ok_cdnc, aerosol_couple, chemistry_couple
     120    INTEGER :: flag_aerosol
     121    INTEGER :: flag_aerosol_strat
     122    INTEGER :: flag_volc_surfstrat
     123    LOGICAL :: flag_aer_feedback
     124    LOGICAL :: flag_bc_internal_mixture
     125    REAL :: bl95_b0, bl95_b1
     126    INTEGER :: read_climoz                        !--- Read ozone climatology
     127    REAL :: alp_offset
     128    LOGICAL :: filtre_oro = .FALSE.
     129
     130    INCLUDE "compbl.h"
     131    INCLUDE "alpale.h"
     132
     133    deg2rad = pi / 180.0
     134    iml = assert_eq(SIZE(masque, 1), SIZE(phis, 1), TRIM(modname) // " iml")
     135    jml = assert_eq(SIZE(masque, 2), SIZE(phis, 2), TRIM(modname) // " jml")
     136
     137    ! Physics configuration
     138    !*******************************************************************************
     139    CALL conf_phys(ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, &
     140            callstats, &
     141            solarlong0, seuil_inversion, &
     142            fact_cldcon, facttemps, ok_newmicro, iflag_radia, &
     143            iflag_cldcon, &
     144            ratqsbas, ratqshaut, tau_ratqs, &
     145            ok_ade, ok_aie, ok_alw, ok_cdnc, ok_volcan, flag_volc_surfstrat, &
     146            aerosol_couple, chemistry_couple, flag_aerosol, flag_aerosol_strat, &
     147            flag_aer_feedback, flag_bc_internal_mixture, bl95_b0, bl95_b1, &
     148            read_climoz, alp_offset)
     149    CALL phys_state_var_init(read_climoz)
     150
     151    !--- Initial atmospheric CO2 conc. from .def file
     152    co2_ppm0 = co2_ppm
     153
     154    ! Compute ground geopotential, sub-cells quantities and possibly the mask.
     155    !*******************************************************************************
     156    read_mask = ANY(masque/=-99999.); masque_tmp = masque
     157    CALL start_init_orog(rlonv, rlatu, phis, masque_tmp)
     158
     159    CALL getin('filtre_oro', filtre_oro)
     160    IF (filtre_oro) CALL filtreoro(size(phis, 1), size(phis, 2), phis, masque_tmp, rlatu)
     161
     162    WRITE(fmt, "(i4,'i1)')")iml ; fmt = '(' // ADJUSTL(fmt)
     163    IF(.NOT.read_mask) THEN                       !--- Keep mask form orography
     164      masque = masque_tmp
     165      IF(prt_level>=1) THEN
     166        WRITE(lunout, *)'BUILT MASK :'
     167        WRITE(lunout, fmt) NINT(masque)
     168      END IF
     169      WHERE(masque(:, :)<EPSFRA) masque(:, :) = 0.
     170      WHERE(1. - masque(:, :)<EPSFRA) masque(:, :) = 1.
    169171    END IF
    170     WHERE(   masque(:,:)<EPSFRA) masque(:,:)=0.
    171     WHERE(1.-masque(:,:)<EPSFRA) masque(:,:)=1.
    172   END IF
    173   CALL gr_dyn_fi(1,iml,jml,klon,masque,zmasq) !--- Land mask to physical grid
    174 
    175 ! Compute tsol and qsol on physical grid, knowing phis on 2D grid.
    176 !*******************************************************************************
    177   CALL start_init_phys(rlonu, rlatv, phis)
    178 
    179 ! Some initializations.
    180 !*******************************************************************************
    181   sn    (:) = 0.0                               !--- Snow
    182   radsol(:) = 0.0                               !--- Net radiation at ground
    183   rugmer(:) = 0.001                             !--- Ocean rugosity
    184   !--- Ozone (zonal or 3D) interpolation in space and time (if 2nd arg is TRUE)
    185   IF(read_climoz>=1) CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
    186 
    187 ! Sub-surfaces initialization.
    188 !*******************************************************************************
    189   CALL start_init_subsurf(read_mask)
    190 
    191 ! Write physical initial state
    192 !*******************************************************************************
    193   WRITE(lunout,*)'phystep ',dtvr,iphysiq,nbapp_rad
    194   phystep = dtvr * FLOAT(iphysiq)
    195   radpas  = NINT (86400./phystep/ FLOAT(nbapp_rad) )
    196   WRITE(lunout,*)'phystep =', phystep, radpas
    197 
    198 ! Init: ftsol, snsrf, qsurf, tsoil, rain_fall, snow_fall, solsw, sollw, z0
    199 !*******************************************************************************
    200   DO i=1,nbsrf; ftsol(:,i) = tsol; END DO
    201   DO i=1,nbsrf; snsrf(:,i) = sn;   END DO
    202   falb_dir(:, :, is_ter) = 0.08
    203   falb_dir(:, :, is_lic) = 0.6
    204   falb_dir(:, :, is_oce) = 0.5
    205   falb_dir(:, :, is_sic) = 0.6
    206 
    207 !ym warning missing init for falb_dif => set to 0
    208   falb_dif(:,:,:)=0
    209 
    210   u10m(:,:)=0 
    211   v10m(:,:)=0 
    212   treedrg(:,:,:)=0
    213 
    214   fevap(:,:) = 0.
    215   qsurf = 0.
    216   DO i=1,nbsrf; DO j=1,nsoilmx; tsoil(:,j,i) = tsol; END DO; END DO
    217   rain_fall  = 0.
    218   snow_fall = 0.
    219   solsw      = 165.
    220   solswfdiff = 1.
    221   sollw      = -53.
    222 !ym warning missing init for sollwdown => set to 0
    223   sollwdown  = 0.
    224   t_ancien   = 273.15
    225   q_ancien   = 0.
    226   ql_ancien = 0.
    227   qs_ancien = 0.
    228   prlw_ancien = 0.
    229   prsw_ancien = 0.
    230   prw_ancien = 0.
    231   agesno    = 0.
    232  
    233   u_ancien = 0.
    234   v_ancien = 0.
    235   wake_delta_pbl_TKE(:,:,:)=0
    236   wake_dens(:)=0
    237   awake_dens = 0.
    238   cv_gen = 0.
    239   ale_bl = 0.
    240   ale_bl_trig =0.
    241   alp_bl=0.
    242   ale_wake=0.
    243   ale_bl_stat=0.
    244  
    245   z0m(:,:)=0 ! ym missing 5th subsurface initialization
    246  
    247   z0m(:,is_oce) = rugmer(:)
    248   z0m(:,is_ter) = 0.01 !MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
    249   z0m(:,is_lic) = 0.001 !MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
    250   z0m(:,is_sic) = 0.001
    251   z0h(:,:)=z0m(:,:)
    252 
    253   fder    = 0.0
    254   clwcon = 0.0
    255   rnebcon = 0.0
    256   ratqs  = 0.0
    257   run_off_lic_0 = 0.0
    258   rugoro = 0.0
    259 
    260 ! Before phyredem calling, surface modules and values to be saved in startphy.nc
    261 ! are initialized
    262 !*******************************************************************************
    263   dummy            = 1.0
    264   pbl_tke(:,:,:)   = 1.e-8
    265   zmax0(:)         = 40.
    266   f0(:)            = 1.e-5
    267   sig1(:,:)        = 0.
    268   w01(:,:)        = 0.
    269   wake_deltat(:,:) = 0.
    270   wake_deltaq(:,:) = 0.
    271   wake_s(:)        = 0.
    272   wake_cstar(:)    = 0.
    273   wake_fip(:)      = 0.
    274   wake_pe          = 0.
    275   fm_therm        = 0.
    276   entr_therm      = 0.
    277   detr_therm      = 0.
    278   awake_s = 0.
    279 
    280   CALL fonte_neige_init(run_off_lic_0)
    281   CALL pbl_surface_init( fder, snsrf, qsurf, tsoil )
    282 
    283   IF (iflag_pbl>1 .AND. iflag_wake>=1  .AND. iflag_pbl_split >=1) THEN
    284      delta_tsurf = 0.
    285      beta_aridity = 0.
    286   end IF
    287 
    288   ratqs_inter_ = 0.002
    289   rneb_ancien = 0.
    290   CALL phyredem( "startphy.nc" )
    291 
    292 !  WRITE(lunout,*)'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0'
    293 !  WRITE(lunout,*)'entree histclo'
    294   CALL histclo()
    295 
    296 END SUBROUTINE etat0phys_netcdf
    297 
    298 !-------------------------------------------------------------------------------
    299 
    300 
    301 !-------------------------------------------------------------------------------
    302 
    303 SUBROUTINE start_init_orog(lon_in,lat_in,phis,masque)
    304 
    305 !===============================================================================
    306 ! Comment:
    307 !   This routine launch grid_noro, which computes parameters for SSO scheme as
    308 !   described in LOTT & MILLER (1997) and LOTT(1999).
    309 !   In case the file oroparam is present and the key read_orop is activated,
    310 !   grid_noro is bypassed and sub-cell parameters are read from the file.
    311 !===============================================================================
    312   USE grid_noro_m, ONLY: grid_noro, read_noro
    313   USE logic_mod,   ONLY: read_orop
    314   IMPLICIT NONE
    315 !-------------------------------------------------------------------------------
    316 ! Arguments:
    317   REAL,    INTENT(IN)    :: lon_in(:), lat_in(:)   ! dim (iml) (jml)
    318   REAL,    INTENT(INOUT) :: phis(:,:), masque(:,:) ! dim (iml,jml)
    319 !-------------------------------------------------------------------------------
    320 ! Local variables:
    321   CHARACTER(LEN=256) :: modname
    322   INTEGER            :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1)
    323   INTEGER            :: ierr
    324   REAL               :: lev(1), date, dt
    325   REAL, ALLOCATABLE  :: lon_rad(:), lon_ini(:), lon_rel(:,:), relief_hi(:,:)
    326   REAL, ALLOCATABLE  :: lat_rad(:), lat_ini(:), lat_rel(:,:), tmp_var  (:,:)
    327   REAL, ALLOCATABLE  :: zmea0(:,:), zstd0(:,:), zsig0(:,:)
    328   REAL, ALLOCATABLE  :: zgam0(:,:), zthe0(:,:), zpic0(:,:), zval0(:,:)
    329 !-------------------------------------------------------------------------------
    330   modname="start_init_orog"
    331   iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml")
    332   jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml")
    333 
    334 !--- HIGH RESOLUTION OROGRAPHY
    335   CALL flininfo(orofname, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
    336 
    337   ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel))
    338   CALL flinopen(orofname, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,&
    339                 lev, ttm_tmp, itau, date, dt, fid)
    340   ALLOCATE(relief_hi(iml_rel,jml_rel))
    341   CALL flinget(fid, orogvar, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1,1, relief_hi)
    342   CALL flinclo(fid)
    343 
    344 !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
    345   ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel))
    346   lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad
    347   lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad
    348 
    349 !--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
    350   ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel))
    351   CALL conf_dat2d(orogvar, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi,.FALSE.)
    352   DEALLOCATE(lon_ini,lat_ini)
    353 
    354 !--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro
    355   WRITE(lunout,*)
    356   WRITE(lunout,*)'*** Compute parameters needed for gravity wave drag code ***'
    357 
    358 !--- ALLOCATIONS OF SUB-CELL SCALES QUANTITIES
    359   ALLOCATE(zmea0(iml,jml),zstd0(iml,jml)) !--- Mean orography and std deviation
    360   ALLOCATE(zsig0(iml,jml),zgam0(iml,jml)) !--- Slope and nisotropy
    361   zsig0(:,:)=0   !ym uninitialized variable
    362   zgam0(:,:)=0   !ym uninitialized variable
    363   ALLOCATE(zthe0(iml,jml))                !--- Highest slope orientation
    364   zthe0(:,:)=0   !ym uninitialized variable
    365   ALLOCATE(zpic0(iml,jml),zval0(iml,jml)) !--- Peaks and valley heights
    366 
    367 !--- READ SUB-CELL SCALES PARAMETERS FROM A FILE (AT RIGHT RESOLUTION)
    368   OPEN(UNIT=66,FILE=oroparam,STATUS='OLD',IOSTAT=ierr)
    369   IF(ierr==0.AND.read_orop) THEN
    370     CLOSE(UNIT=66)
    371     CALL read_noro(lon_in,lat_in,oroparam,                                     &
    372                    phis,zmea0,zstd0,zsig0,zgam0,zthe0,zpic0,zval0,masque)
    373   ELSE
    374 !--- CALL OROGRAPHY MODULE TO COMPUTE FIELDS
    375     CALL grid_noro(lon_rad,lat_rad,relief_hi,lon_in,lat_in,                    &
    376                    phis,zmea0,zstd0,zsig0,zgam0,zthe0,zpic0,zval0,masque)
    377   END IF
    378   phis = phis * 9.81
    379   phis(iml,:) = phis(1,:)
    380   DEALLOCATE(relief_hi,lon_rad,lat_rad)
    381 
    382 !--- PUT QUANTITIES TO PHYSICAL GRID
    383   CALL gr_dyn_fi(1,iml,jml,klon,zmea0,zmea); DEALLOCATE(zmea0)
    384   CALL gr_dyn_fi(1,iml,jml,klon,zstd0,zstd); DEALLOCATE(zstd0)
    385   CALL gr_dyn_fi(1,iml,jml,klon,zsig0,zsig); DEALLOCATE(zsig0)
    386   CALL gr_dyn_fi(1,iml,jml,klon,zgam0,zgam); DEALLOCATE(zgam0)
    387   CALL gr_dyn_fi(1,iml,jml,klon,zthe0,zthe); DEALLOCATE(zthe0)
    388   CALL gr_dyn_fi(1,iml,jml,klon,zpic0,zpic); DEALLOCATE(zpic0)
    389   CALL gr_dyn_fi(1,iml,jml,klon,zval0,zval); DEALLOCATE(zval0)
    390 
    391 
    392 END SUBROUTINE start_init_orog
    393 
    394 !-------------------------------------------------------------------------------
    395 
    396 
    397 !-------------------------------------------------------------------------------
    398 
    399 SUBROUTINE start_init_phys(lon_in,lat_in,phis)
    400 
    401 !===============================================================================
    402 ! Purpose:   Compute tsol and qsol, knowing phis.
    403 !===============================================================================
    404   IMPLICIT NONE
    405 !-------------------------------------------------------------------------------
    406 ! Arguments:
    407   REAL,    INTENT(IN) :: lon_in(:), lat_in(:)       ! dim (iml) (jml2)
    408   REAL,    INTENT(IN) :: phis(:,:)                   ! dim (iml,jml)
    409 !-------------------------------------------------------------------------------
    410 ! Local variables:
    411   CHARACTER(LEN=256) :: modname
    412   REAL              :: date, dt
    413   INTEGER            :: iml, jml, jml2, itau(1)
    414   REAL, ALLOCATABLE  :: lon_rad(:), lon_ini(:), var_ana(:,:)
    415   REAL, ALLOCATABLE :: lat_rad(:), lat_ini(:)
    416   REAL, ALLOCATABLE  :: ts(:,:), qs(:,:)
    417 !-------------------------------------------------------------------------------
    418   modname="start_init_phys"
    419   iml=assert_eq(SIZE(lon_in),SIZE(phis,1),TRIM(modname)//" iml")
    420   jml=SIZE(phis,2); jml2=SIZE(lat_in)
    421 
    422   WRITE(lunout,*)'Opening the surface analysis'
    423   CALL flininfo(phyfname, iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys)
    424   WRITE(lunout,*) 'Values read: ', iml_phys, jml_phys, llm_phys, ttm_phys
    425 
    426   ALLOCATE(lat_phys(iml_phys,jml_phys),lon_phys(iml_phys,jml_phys))
    427   ALLOCATE(levphys_ini(llm_phys))
    428   CALL flinopen(phyfname, .FALSE., iml_phys, jml_phys, llm_phys,              &
    429                 lon_phys,lat_phys,levphys_ini,ttm_phys,itau,date,dt,fid_phys)
    430 
    431 !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
    432   ALLOCATE(lon_ini(iml_phys),lat_ini(jml_phys))
    433   lon_ini(:)=lon_phys(:,1); IF(MAXVAL(lon_phys)>pi) lon_ini=lon_ini*deg2rad
    434   lat_ini(:)=lat_phys(1,:); IF(MAXVAL(lat_phys)>pi) lat_ini=lat_ini*deg2rad
    435 
    436   ALLOCATE(var_ana(iml_phys,jml_phys),lon_rad(iml_phys),lat_rad(jml_phys))
    437   CALL get_var_phys(tsrfvar,ts)                   !--- SURFACE TEMPERATURE
    438   CALL get_var_phys(qsolvar,qs)                   !--- SOIL MOISTURE
    439   CALL flinclo(fid_phys)
    440   DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
    441 
    442 !--- TSOL AND QSOL ON PHYSICAL GRID
    443   ALLOCATE(tsol(klon))
    444   CALL gr_dyn_fi(1,iml,jml,klon,ts,tsol)
    445   CALL gr_dyn_fi(1,iml,jml,klon,qs,qsol)
    446   DEALLOCATE(ts,qs)
    447 
    448 CONTAINS
    449 
    450 !-------------------------------------------------------------------------------
    451 
    452 SUBROUTINE get_var_phys(title,field)
    453 
    454 !-------------------------------------------------------------------------------
    455   IMPLICIT NONE
    456 !-------------------------------------------------------------------------------
    457 ! Arguments:
    458   CHARACTER(LEN=*),  INTENT(IN)    :: title
    459   REAL, ALLOCATABLE, INTENT(INOUT) :: field(:,:)
    460 !-------------------------------------------------------------------------------
    461 ! Local variables:
    462   INTEGER :: tllm
    463 !-------------------------------------------------------------------------------
    464   SELECT CASE(title)
    465     CASE(psrfvar);         tllm=0
    466     CASE(tsrfvar,qsolvar); tllm=llm_phys
    467   END SELECT
    468   IF(ALLOCATED(field)) RETURN
    469   ALLOCATE(field(iml,jml)); field(:,:)=0.
    470   CALL flinget(fid_phys,title,iml_phys,jml_phys,tllm,ttm_phys,1,1,var_ana)
    471   CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, .TRUE.)
    472   CALL interp_startvar(title, .TRUE., lon_rad, lat_rad, var_ana,              &
    473                                       lon_in,  lat_in, field)
    474 
    475 END SUBROUTINE get_var_phys
    476 
    477 !-------------------------------------------------------------------------------
    478 
    479 END SUBROUTINE start_init_phys
    480 
    481 !-------------------------------------------------------------------------------
    482 
    483 
    484 !-------------------------------------------------------------------------------
    485 
    486 SUBROUTINE interp_startvar(nam,ibeg,lon,lat,vari,lon2,lat2,varo)
    487 
    488 !-------------------------------------------------------------------------------
    489   USE inter_barxy_m, ONLY: inter_barxy
    490   IMPLICIT NONE
    491 !-------------------------------------------------------------------------------
    492 ! Arguments:
    493   CHARACTER(LEN=*), INTENT(IN) :: nam
    494   LOGICAL,          INTENT(IN) :: ibeg
    495   REAL,             INTENT(IN) :: lon(:), lat(:)   ! dim (ii) (jj)
    496   REAL,             INTENT(IN)  :: vari(:,:)        ! dim (ii,jj)
    497   REAL,             INTENT(IN) :: lon2(:), lat2(:) ! dim (i1) (j2)
    498   REAL,             INTENT(OUT) :: varo(:,:)        ! dim (i1) (j1)
    499 !-------------------------------------------------------------------------------
    500 ! Local variables:
    501   CHARACTER(LEN=256) :: modname
    502   INTEGER            :: ii, jj, i1, j1, j2
    503   REAL, ALLOCATABLE  :: vtmp(:,:)
    504 !-------------------------------------------------------------------------------
    505   modname="interp_startvar"
    506   ii=assert_eq(SIZE(lon), SIZE(vari,1),TRIM(modname)//" ii")
    507   jj=assert_eq(SIZE(lat), SIZE(vari,2),TRIM(modname)//" jj")
    508   i1=assert_eq(SIZE(lon2),SIZE(varo,1),TRIM(modname)//" i1")
    509   j1=SIZE(varo,2); j2=SIZE(lat2)
    510   ALLOCATE(vtmp(i1-1,j1))
    511   IF(ibeg.AND.prt_level>1) THEN
    512     WRITE(lunout,*)"--------------------------------------------------------"
    513     WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$"
    514     WRITE(lunout,*)"--------------------------------------------------------"
    515   END IF
    516   CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp)
    517   CALL gr_int_dyn(vtmp, varo, i1-1, j1)
    518 
    519 END SUBROUTINE interp_startvar
    520 
    521 !-------------------------------------------------------------------------------
    522 
    523 !*******************************************************************************
    524 
    525 SUBROUTINE filtreoro(imp1,jmp1,phis,masque,rlatu)
    526 
    527 IMPLICIT NONE
    528 
    529   INTEGER imp1,jmp1
    530   REAL, DIMENSION(imp1,jmp1) :: phis,masque
    531   REAL, DIMENSION(jmp1) :: rlatu
    532   REAL, DIMENSION(imp1) :: wwf
    533   REAL, DIMENSION(imp1,jmp1) :: phiso
    534   INTEGER :: ifiltre,ifi,ii,i,j
    535   REAL :: coslat0,ssz
    536 
    537   coslat0=0.5
    538   phiso=phis
    539   do j=2,jmp1-1
    540      PRINT*,'avant if ',cos(rlatu(j)),coslat0
    541      IF (cos(rlatu(j))<coslat0) THEN
    542          ! nb de pts affectes par le filtrage de part et d'autre du pt
    543          ifiltre=(coslat0/cos(rlatu(j))-1.)/2.
    544          wwf=0.
    545          do i=1,ifiltre
    546             wwf(i)=1.
    547          enddo
    548          wwf(ifiltre+1)=(coslat0/cos(rlatu(j))-1.)/2.-ifiltre
    549          do i=1,imp1-1
    550             IF (masque(i,j)>0.9) THEN
    551                ssz=phis(i,j)
    552                do ifi=1,ifiltre+1
    553                   ii=i+ifi
    554                   IF (ii>imp1-1) ii=ii-imp1+1
    555                   ssz=ssz+wwf(ifi)*phis(ii,j)
    556                   ii=i-ifi
    557                   IF (ii<1) ii=ii+imp1-1
    558                   ssz=ssz+wwf(ifi)*phis(ii,j)
    559                enddo
    560                phis(i,j)=ssz*cos(rlatu(j))/coslat0
    561             endif
    562          enddo
    563          PRINT*,'j=',j,coslat0/cos(rlatu(j)), (1.+2.*sum(wwf))*cos(rlatu(j))/coslat0
    564      endif
    565   enddo
    566   CALL dump2d(imp1,jmp1,phis,'phis ')
    567   CALL dump2d(imp1,jmp1,masque,'masque ')
    568   CALL dump2d(imp1,jmp1,phis-phiso,'dphis ')
    569 
    570 END SUBROUTINE filtreoro
     172    CALL gr_dyn_fi(1, iml, jml, klon, masque, zmasq) !--- Land mask to physical grid
     173
     174    ! Compute tsol and qsol on physical grid, knowing phis on 2D grid.
     175    !*******************************************************************************
     176    CALL start_init_phys(rlonu, rlatv, phis)
     177
     178    ! Some initializations.
     179    !*******************************************************************************
     180    sn    (:) = 0.0                               !--- Snow
     181    radsol(:) = 0.0                               !--- Net radiation at ground
     182    rugmer(:) = 0.001                             !--- Ocean rugosity
     183    !--- Ozone (zonal or 3D) interpolation in space and time (if 2nd arg is TRUE)
     184    IF(read_climoz>=1) CALL regr_horiz_time_climoz(read_climoz, ok_daily_climoz)
     185
     186    ! Sub-surfaces initialization.
     187    !*******************************************************************************
     188    CALL start_init_subsurf(read_mask)
     189
     190    ! Write physical initial state
     191    !*******************************************************************************
     192    WRITE(lunout, *)'phystep ', dtvr, iphysiq, nbapp_rad
     193    phystep = dtvr * FLOAT(iphysiq)
     194    radpas = NINT (86400. / phystep / FLOAT(nbapp_rad))
     195    WRITE(lunout, *)'phystep =', phystep, radpas
     196
     197    ! Init: ftsol, snsrf, qsurf, tsoil, rain_fall, snow_fall, solsw, sollw, z0
     198    !*******************************************************************************
     199    DO i = 1, nbsrf; ftsol(:, i) = tsol;
     200    END DO
     201    DO i = 1, nbsrf; snsrf(:, i) = sn;
     202    END DO
     203    falb_dir(:, :, is_ter) = 0.08
     204    falb_dir(:, :, is_lic) = 0.6
     205    falb_dir(:, :, is_oce) = 0.5
     206    falb_dir(:, :, is_sic) = 0.6
     207
     208    !ym warning missing init for falb_dif => set to 0
     209    falb_dif(:, :, :) = 0
     210
     211    u10m(:, :) = 0
     212    v10m(:, :) = 0
     213    treedrg(:, :, :) = 0
     214
     215    fevap(:, :) = 0.
     216    qsurf = 0.
     217    DO i = 1, nbsrf; DO j = 1, nsoilmx; tsoil(:, j, i) = tsol;
     218    END DO;
     219    END DO
     220    rain_fall = 0.
     221    snow_fall = 0.
     222    solsw = 165.
     223    solswfdiff = 1.
     224    sollw = -53.
     225    !ym warning missing init for sollwdown => set to 0
     226    sollwdown = 0.
     227    t_ancien = 273.15
     228    q_ancien = 0.
     229    ql_ancien = 0.
     230    qs_ancien = 0.
     231    prlw_ancien = 0.
     232    prsw_ancien = 0.
     233    prw_ancien = 0.
     234    agesno = 0.
     235
     236    u_ancien = 0.
     237    v_ancien = 0.
     238    wake_delta_pbl_TKE(:, :, :) = 0
     239    wake_dens(:) = 0
     240    awake_dens = 0.
     241    cv_gen = 0.
     242    ale_bl = 0.
     243    ale_bl_trig = 0.
     244    alp_bl = 0.
     245    ale_wake = 0.
     246    ale_bl_stat = 0.
     247
     248    z0m(:, :) = 0 ! ym missing 5th subsurface initialization
     249
     250    z0m(:, is_oce) = rugmer(:)
     251    z0m(:, is_ter) = 0.01 !MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
     252    z0m(:, is_lic) = 0.001 !MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
     253    z0m(:, is_sic) = 0.001
     254    z0h(:, :) = z0m(:, :)
     255
     256    fder = 0.0
     257    clwcon = 0.0
     258    rnebcon = 0.0
     259    ratqs = 0.0
     260    run_off_lic_0 = 0.0
     261    rugoro = 0.0
     262
     263    ! Before phyredem calling, surface modules and values to be saved in startphy.nc
     264    ! are initialized
     265    !*******************************************************************************
     266    dummy = 1.0
     267    pbl_tke(:, :, :) = 1.e-8
     268    zmax0(:) = 40.
     269    f0(:) = 1.e-5
     270    sig1(:, :) = 0.
     271    w01(:, :) = 0.
     272    wake_deltat(:, :) = 0.
     273    wake_deltaq(:, :) = 0.
     274    wake_s(:) = 0.
     275    wake_cstar(:) = 0.
     276    wake_fip(:) = 0.
     277    wake_pe = 0.
     278    fm_therm = 0.
     279    entr_therm = 0.
     280    detr_therm = 0.
     281    awake_s = 0.
     282
     283    CALL fonte_neige_init(run_off_lic_0)
     284    CALL pbl_surface_init(fder, snsrf, qsurf, tsoil)
     285
     286    IF (iflag_pbl>1 .AND. iflag_wake>=1  .AND. iflag_pbl_split >=1) THEN
     287      delta_tsurf = 0.
     288      beta_aridity = 0.
     289    end IF
     290
     291    ratqs_inter_ = 0.002
     292    rneb_ancien = 0.
     293    CALL phyredem("startphy.nc")
     294
     295    !  WRITE(lunout,*)'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0'
     296    !  WRITE(lunout,*)'entree histclo'
     297    CALL histclo()
     298
     299  END SUBROUTINE etat0phys_netcdf
     300
     301  !-------------------------------------------------------------------------------
     302
     303
     304  !-------------------------------------------------------------------------------
     305
     306  SUBROUTINE start_init_orog(lon_in, lat_in, phis, masque)
     307
     308    !===============================================================================
     309    ! Comment:
     310    !   This routine launch grid_noro, which computes parameters for SSO scheme as
     311    !   described in LOTT & MILLER (1997) and LOTT(1999).
     312    !   In case the file oroparam is present and the key read_orop is activated,
     313    !   grid_noro is bypassed and sub-cell parameters are read from the file.
     314    !===============================================================================
     315    USE grid_noro_m, ONLY: grid_noro, read_noro
     316    USE logic_mod, ONLY: read_orop
     317    IMPLICIT NONE
     318    !-------------------------------------------------------------------------------
     319    ! Arguments:
     320    REAL, INTENT(IN) :: lon_in(:), lat_in(:)   ! dim (iml) (jml)
     321    REAL, INTENT(INOUT) :: phis(:, :), masque(:, :) ! dim (iml,jml)
     322    !-------------------------------------------------------------------------------
     323    ! Local variables:
     324    CHARACTER(LEN = 256) :: modname
     325    INTEGER :: fid, llm_tmp, ttm_tmp, iml, jml, iml_rel, jml_rel, itau(1)
     326    INTEGER :: ierr
     327    REAL :: lev(1), date, dt
     328    REAL, ALLOCATABLE :: lon_rad(:), lon_ini(:), lon_rel(:, :), relief_hi(:, :)
     329    REAL, ALLOCATABLE :: lat_rad(:), lat_ini(:), lat_rel(:, :), tmp_var  (:, :)
     330    REAL, ALLOCATABLE :: zmea0(:, :), zstd0(:, :), zsig0(:, :)
     331    REAL, ALLOCATABLE :: zgam0(:, :), zthe0(:, :), zpic0(:, :), zval0(:, :)
     332    !-------------------------------------------------------------------------------
     333    modname = "start_init_orog"
     334    iml = assert_eq(SIZE(lon_in), SIZE(phis, 1), SIZE(masque, 1), TRIM(modname) // " iml")
     335    jml = assert_eq(SIZE(lat_in), SIZE(phis, 2), SIZE(masque, 2), TRIM(modname) // " jml")
     336
     337    !--- HIGH RESOLUTION OROGRAPHY
     338    CALL flininfo(orofname, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
     339
     340    ALLOCATE(lat_rel(iml_rel, jml_rel), lon_rel(iml_rel, jml_rel))
     341    CALL flinopen(orofname, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel, &
     342            lev, ttm_tmp, itau, date, dt, fid)
     343    ALLOCATE(relief_hi(iml_rel, jml_rel))
     344    CALL flinget(fid, orogvar, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi)
     345    CALL flinclo(fid)
     346
     347    !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
     348    ALLOCATE(lon_ini(iml_rel), lat_ini(jml_rel))
     349    lon_ini(:) = lon_rel(:, 1); IF(MAXVAL(lon_rel)>pi) lon_ini = lon_ini * deg2rad
     350    lat_ini(:) = lat_rel(1, :); IF(MAXVAL(lat_rel)>pi) lat_ini = lat_ini * deg2rad
     351
     352    !--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
     353    ALLOCATE(lon_rad(iml_rel), lat_rad(jml_rel))
     354    CALL conf_dat2d(orogvar, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi, .FALSE.)
     355    DEALLOCATE(lon_ini, lat_ini)
     356
     357    !--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro
     358    WRITE(lunout, *)
     359    WRITE(lunout, *)'*** Compute parameters needed for gravity wave drag code ***'
     360
     361    !--- ALLOCATIONS OF SUB-CELL SCALES QUANTITIES
     362    ALLOCATE(zmea0(iml, jml), zstd0(iml, jml)) !--- Mean orography and std deviation
     363    ALLOCATE(zsig0(iml, jml), zgam0(iml, jml)) !--- Slope and nisotropy
     364    zsig0(:, :) = 0   !ym uninitialized variable
     365    zgam0(:, :) = 0   !ym uninitialized variable
     366    ALLOCATE(zthe0(iml, jml))                !--- Highest slope orientation
     367    zthe0(:, :) = 0   !ym uninitialized variable
     368    ALLOCATE(zpic0(iml, jml), zval0(iml, jml)) !--- Peaks and valley heights
     369
     370    !--- READ SUB-CELL SCALES PARAMETERS FROM A FILE (AT RIGHT RESOLUTION)
     371    OPEN(UNIT = 66, FILE = oroparam, STATUS = 'OLD', IOSTAT = ierr)
     372    IF(ierr==0.AND.read_orop) THEN
     373      CLOSE(UNIT = 66)
     374      CALL read_noro(lon_in, lat_in, oroparam, &
     375              phis, zmea0, zstd0, zsig0, zgam0, zthe0, zpic0, zval0, masque)
     376    ELSE
     377      !--- CALL OROGRAPHY MODULE TO COMPUTE FIELDS
     378      CALL grid_noro(lon_rad, lat_rad, relief_hi, lon_in, lat_in, &
     379              phis, zmea0, zstd0, zsig0, zgam0, zthe0, zpic0, zval0, masque)
     380    END IF
     381    phis = phis * 9.81
     382    phis(iml, :) = phis(1, :)
     383    DEALLOCATE(relief_hi, lon_rad, lat_rad)
     384
     385    !--- PUT QUANTITIES TO PHYSICAL GRID
     386    CALL gr_dyn_fi(1, iml, jml, klon, zmea0, zmea); DEALLOCATE(zmea0)
     387    CALL gr_dyn_fi(1, iml, jml, klon, zstd0, zstd); DEALLOCATE(zstd0)
     388    CALL gr_dyn_fi(1, iml, jml, klon, zsig0, zsig); DEALLOCATE(zsig0)
     389    CALL gr_dyn_fi(1, iml, jml, klon, zgam0, zgam); DEALLOCATE(zgam0)
     390    CALL gr_dyn_fi(1, iml, jml, klon, zthe0, zthe); DEALLOCATE(zthe0)
     391    CALL gr_dyn_fi(1, iml, jml, klon, zpic0, zpic); DEALLOCATE(zpic0)
     392    CALL gr_dyn_fi(1, iml, jml, klon, zval0, zval); DEALLOCATE(zval0)
     393
     394  END SUBROUTINE start_init_orog
     395
     396  !-------------------------------------------------------------------------------
     397
     398
     399  !-------------------------------------------------------------------------------
     400
     401  SUBROUTINE start_init_phys(lon_in, lat_in, phis)
     402
     403    !===============================================================================
     404    ! Purpose:   Compute tsol and qsol, knowing phis.
     405    !===============================================================================
     406    IMPLICIT NONE
     407    !-------------------------------------------------------------------------------
     408    ! Arguments:
     409    REAL, INTENT(IN) :: lon_in(:), lat_in(:)       ! dim (iml) (jml2)
     410    REAL, INTENT(IN) :: phis(:, :)                   ! dim (iml,jml)
     411    !-------------------------------------------------------------------------------
     412    ! Local variables:
     413    CHARACTER(LEN = 256) :: modname
     414    REAL :: date, dt
     415    INTEGER :: iml, jml, jml2, itau(1)
     416    REAL, ALLOCATABLE :: lon_rad(:), lon_ini(:), var_ana(:, :)
     417    REAL, ALLOCATABLE :: lat_rad(:), lat_ini(:)
     418    REAL, ALLOCATABLE :: ts(:, :), qs(:, :)
     419    !-------------------------------------------------------------------------------
     420    modname = "start_init_phys"
     421    iml = assert_eq(SIZE(lon_in), SIZE(phis, 1), TRIM(modname) // " iml")
     422    jml = SIZE(phis, 2); jml2 = SIZE(lat_in)
     423
     424    WRITE(lunout, *)'Opening the surface analysis'
     425    CALL flininfo(phyfname, iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys)
     426    WRITE(lunout, *) 'Values read: ', iml_phys, jml_phys, llm_phys, ttm_phys
     427
     428    ALLOCATE(lat_phys(iml_phys, jml_phys), lon_phys(iml_phys, jml_phys))
     429    ALLOCATE(levphys_ini(llm_phys))
     430    CALL flinopen(phyfname, .FALSE., iml_phys, jml_phys, llm_phys, &
     431            lon_phys, lat_phys, levphys_ini, ttm_phys, itau, date, dt, fid_phys)
     432
     433    !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
     434    ALLOCATE(lon_ini(iml_phys), lat_ini(jml_phys))
     435    lon_ini(:) = lon_phys(:, 1); IF(MAXVAL(lon_phys)>pi) lon_ini = lon_ini * deg2rad
     436    lat_ini(:) = lat_phys(1, :); IF(MAXVAL(lat_phys)>pi) lat_ini = lat_ini * deg2rad
     437
     438    ALLOCATE(var_ana(iml_phys, jml_phys), lon_rad(iml_phys), lat_rad(jml_phys))
     439    CALL get_var_phys(tsrfvar, ts)                   !--- SURFACE TEMPERATURE
     440    CALL get_var_phys(qsolvar, qs)                   !--- SOIL MOISTURE
     441    CALL flinclo(fid_phys)
     442    DEALLOCATE(var_ana, lon_rad, lat_rad, lon_ini, lat_ini)
     443
     444    !--- TSOL AND QSOL ON PHYSICAL GRID
     445    ALLOCATE(tsol(klon))
     446    CALL gr_dyn_fi(1, iml, jml, klon, ts, tsol)
     447    CALL gr_dyn_fi(1, iml, jml, klon, qs, qsol)
     448    DEALLOCATE(ts, qs)
     449
     450  CONTAINS
     451
     452    !-------------------------------------------------------------------------------
     453
     454    SUBROUTINE get_var_phys(title, field)
     455
     456      !-------------------------------------------------------------------------------
     457      IMPLICIT NONE
     458      !-------------------------------------------------------------------------------
     459      ! Arguments:
     460      CHARACTER(LEN = *), INTENT(IN) :: title
     461      REAL, ALLOCATABLE, INTENT(INOUT) :: field(:, :)
     462      !-------------------------------------------------------------------------------
     463      ! Local variables:
     464      INTEGER :: tllm
     465      !-------------------------------------------------------------------------------
     466      SELECT CASE(title)
     467      CASE(psrfvar);         tllm = 0
     468      CASE(tsrfvar, qsolvar); tllm = llm_phys
     469      END SELECT
     470      IF(ALLOCATED(field)) RETURN
     471      ALLOCATE(field(iml, jml)); field(:, :) = 0.
     472      CALL flinget(fid_phys, title, iml_phys, jml_phys, tllm, ttm_phys, 1, 1, var_ana)
     473      CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, .TRUE.)
     474      CALL interp_startvar(title, .TRUE., lon_rad, lat_rad, var_ana, &
     475              lon_in, lat_in, field)
     476
     477    END SUBROUTINE get_var_phys
     478
     479    !-------------------------------------------------------------------------------
     480
     481  END SUBROUTINE start_init_phys
     482
     483  !-------------------------------------------------------------------------------
     484
     485
     486  !-------------------------------------------------------------------------------
     487
     488  SUBROUTINE interp_startvar(nam, ibeg, lon, lat, vari, lon2, lat2, varo)
     489
     490    !-------------------------------------------------------------------------------
     491    USE inter_barxy_m, ONLY: inter_barxy
     492    IMPLICIT NONE
     493    !-------------------------------------------------------------------------------
     494    ! Arguments:
     495    CHARACTER(LEN = *), INTENT(IN) :: nam
     496    LOGICAL, INTENT(IN) :: ibeg
     497    REAL, INTENT(IN) :: lon(:), lat(:)   ! dim (ii) (jj)
     498    REAL, INTENT(IN) :: vari(:, :)        ! dim (ii,jj)
     499    REAL, INTENT(IN) :: lon2(:), lat2(:) ! dim (i1) (j2)
     500    REAL, INTENT(OUT) :: varo(:, :)        ! dim (i1) (j1)
     501    !-------------------------------------------------------------------------------
     502    ! Local variables:
     503    CHARACTER(LEN = 256) :: modname
     504    INTEGER :: ii, jj, i1, j1, j2
     505    REAL, ALLOCATABLE :: vtmp(:, :)
     506    !-------------------------------------------------------------------------------
     507    modname = "interp_startvar"
     508    ii = assert_eq(SIZE(lon), SIZE(vari, 1), TRIM(modname) // " ii")
     509    jj = assert_eq(SIZE(lat), SIZE(vari, 2), TRIM(modname) // " jj")
     510    i1 = assert_eq(SIZE(lon2), SIZE(varo, 1), TRIM(modname) // " i1")
     511    j1 = SIZE(varo, 2); j2 = SIZE(lat2)
     512    ALLOCATE(vtmp(i1 - 1, j1))
     513    IF(ibeg.AND.prt_level>1) THEN
     514      WRITE(lunout, *)"--------------------------------------------------------"
     515      WRITE(lunout, *)"$$$ Interpolation barycentrique pour " // TRIM(nam) // " $$$"
     516      WRITE(lunout, *)"--------------------------------------------------------"
     517    END IF
     518    CALL inter_barxy(lon, lat(:jj - 1), vari, lon2(:i1 - 1), lat2, vtmp)
     519    CALL gr_int_dyn(vtmp, varo, i1 - 1, j1)
     520
     521  END SUBROUTINE interp_startvar
     522
     523  !-------------------------------------------------------------------------------
     524
     525  !*******************************************************************************
     526
     527  SUBROUTINE filtreoro(imp1, jmp1, phis, masque, rlatu)
     528
     529    IMPLICIT NONE
     530
     531    INTEGER imp1, jmp1
     532    REAL, DIMENSION(imp1, jmp1) :: phis, masque
     533    REAL, DIMENSION(jmp1) :: rlatu
     534    REAL, DIMENSION(imp1) :: wwf
     535    REAL, DIMENSION(imp1, jmp1) :: phiso
     536    INTEGER :: ifiltre, ifi, ii, i, j
     537    REAL :: coslat0, ssz
     538
     539    coslat0 = 0.5
     540    phiso = phis
     541    do j = 2, jmp1 - 1
     542      PRINT*, 'avant if ', cos(rlatu(j)), coslat0
     543      IF (cos(rlatu(j))<coslat0) THEN
     544        ! nb de pts affectes par le filtrage de part et d'autre du pt
     545        ifiltre = (coslat0 / cos(rlatu(j)) - 1.) / 2.
     546        wwf = 0.
     547        do i = 1, ifiltre
     548          wwf(i) = 1.
     549        enddo
     550        wwf(ifiltre + 1) = (coslat0 / cos(rlatu(j)) - 1.) / 2. - ifiltre
     551        do i = 1, imp1 - 1
     552          IF (masque(i, j)>0.9) THEN
     553            ssz = phis(i, j)
     554            do ifi = 1, ifiltre + 1
     555              ii = i + ifi
     556              IF (ii>imp1 - 1) ii = ii - imp1 + 1
     557              ssz = ssz + wwf(ifi) * phis(ii, j)
     558              ii = i - ifi
     559              IF (ii<1) ii = ii + imp1 - 1
     560              ssz = ssz + wwf(ifi) * phis(ii, j)
     561            enddo
     562            phis(i, j) = ssz * cos(rlatu(j)) / coslat0
     563          endif
     564        enddo
     565        PRINT*, 'j=', j, coslat0 / cos(rlatu(j)), (1. + 2. * sum(wwf)) * cos(rlatu(j)) / coslat0
     566      endif
     567    enddo
     568    CALL dump2d(imp1, jmp1, phis, 'phis ')
     569    CALL dump2d(imp1, jmp1, masque, 'masque ')
     570    CALL dump2d(imp1, jmp1, phis - phiso, 'dphis ')
     571
     572  END SUBROUTINE filtreoro
    571573
    572574

  • LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90

    r5113 r5118  
    1 
    21! $Id$
    32
    … …  
    65CONTAINS
    76
    8 SUBROUTINE iniphysiq(ii,jj,nlayer, &
    9                      nbp, communicator, &
    10                      punjours, pdayref,ptimestep, &
    11                      rlatudyn,rlatvdyn,rlonudyn,rlonvdyn,airedyn,cudyn,cvdyn, &
    12                      prad,pg,pr,pcpp,iflag_phys)
    13   USE dimphy, ONLY: init_dimphy
    14   USE inigeomphy_mod, ONLY: inigeomphy
    15   USE lmdz_grid_phy, ONLY: nbp_lon,nbp_lat,nbp_lev,klon_glo ! number of atmospheric columns (on full grid)
    16   USE lmdz_phys_para, ONLY: klon_omp ! number of columns (on local omp grid)
    17   USE lmdz_vertical_layers, ONLY: init_vertical_layers
    18   USE infotrac, ONLY: nbtr, type_trac
     7  SUBROUTINE iniphysiq(ii, jj, nlayer, &
     8          nbp, communicator, &
     9          punjours, pdayref, ptimestep, &
     10          rlatudyn, rlatvdyn, rlonudyn, rlonvdyn, airedyn, cudyn, cvdyn, &
     11          prad, pg, pr, pcpp, iflag_phys)
     12    USE dimphy, ONLY: init_dimphy
     13    USE inigeomphy_mod, ONLY: inigeomphy
     14    USE lmdz_grid_phy, ONLY: nbp_lon, nbp_lat, nbp_lev, klon_glo ! number of atmospheric columns (on full grid)
     15    USE lmdz_phys_para, ONLY: klon_omp ! number of columns (on local omp grid)
     16    USE lmdz_vertical_layers, ONLY: init_vertical_layers
     17    USE infotrac, ONLY: nbtr, type_trac
    1918
    2019#ifdef REPROBUS
    … …  
    2625  USE lmdz_phys_omp_data, ONLY: klon_omp
    2726#endif
    28   USE control_mod, ONLY: dayref,anneeref,day_step,nday,offline, iphysiq
    29   USE inifis_mod, ONLY: inifis
    30   USE time_phylmdz_mod, ONLY: init_time
    31   USE temps_mod, ONLY: annee_ref, day_ini, day_ref, start_time, calend, year_len
    32   USE infotrac_phy, ONLY: init_infotrac_phy
    33   USE phystokenc_mod, ONLY: init_phystokenc
    34   USE phyaqua_mod, ONLY: iniaqua
    35   USE comconst_mod, ONLY: omeg, rad
    36   USE indice_sol_mod, ONLY: nbsrf, is_oce, is_sic, is_ter, is_lic
     27    USE control_mod, ONLY: dayref, anneeref, day_step, nday, offline, iphysiq
     28    USE inifis_mod, ONLY: inifis
     29    USE time_phylmdz_mod, ONLY: init_time
     30    USE temps_mod, ONLY: annee_ref, day_ini, day_ref, start_time, calend, year_len
     31    USE infotrac_phy, ONLY: init_infotrac_phy
     32    USE phystokenc_mod, ONLY: init_phystokenc
     33    USE phyaqua_mod, ONLY: iniaqua
     34    USE comconst_mod, ONLY: omeg, rad
     35    USE indice_sol_mod, ONLY: nbsrf, is_oce, is_sic, is_ter, is_lic
    3736#ifdef CPP_PARA
    3837  USE parallel_lmdz, ONLY: mpi_size, mpi_rank
    3938  USE bands, ONLY: distrib_phys
    4039#endif
    41   USE lmdz_phys_omp_data, ONLY: klon_omp
    42   USE lmdz_ioipsl_getin_p, ONLY: getin_p
    43   USE slab_heat_transp_mod, ONLY: ini_slab_transp_geom
    44   IMPLICIT NONE
     40    USE lmdz_phys_omp_data, ONLY: klon_omp
     41    USE lmdz_ioipsl_getin_p, ONLY: getin_p
     42    USE slab_heat_transp_mod, ONLY: ini_slab_transp_geom
     43    USE lmdz_iniprint, ONLY: lunout, prt_level
     44    IMPLICIT NONE
    4545
    46   ! =======================================================================
    47   ! Initialisation of the physical constants and some positional and
    48   ! geometrical arrays for the physics
    49   ! =======================================================================
     46    ! =======================================================================
     47    ! Initialisation of the physical constants and some positional and
     48    ! geometrical arrays for the physics
     49    ! =======================================================================
    5050
    51   include "dimensions.h"
    52   include "paramet.h"
    53   include "iniprint.h"
    54   include "tracstoke.h"
    55   include "comgeom.h"
     51    include "dimensions.h"
     52    include "paramet.h"
     53    include "tracstoke.h"
     54    include "comgeom.h"
    5655
    57   REAL, INTENT (IN) :: prad ! radius of the planet (m)
    58   REAL, INTENT (IN) :: pg ! gravitational acceleration (m/s2)
    59   REAL, INTENT (IN) :: pr ! reduced gas constant R/mu
    60   REAL, INTENT (IN) :: pcpp ! specific heat Cp
    61   REAL, INTENT (IN) :: punjours ! length (in s) of a standard day
    62   INTEGER, INTENT (IN) :: nlayer ! number of atmospheric layers
    63   INTEGER, INTENT (IN) :: ii ! number of atmospheric columns along longitudes
    64   INTEGER, INTENT (IN) :: jj ! number of atompsheric columns along latitudes
    65   INTEGER, INTENT(IN) :: nbp ! number of physics columns for this MPI process
    66   INTEGER, INTENT(IN) :: communicator ! MPI communicator
    67   REAL, INTENT (IN) :: rlatudyn(jj+1) ! latitudes of the physics grid
    68   REAL, INTENT (IN) :: rlatvdyn(jj) ! latitude boundaries of the physics grid
    69   REAL, INTENT (IN) :: rlonvdyn(ii+1) ! longitudes of the physics grid
    70   REAL, INTENT (IN) :: rlonudyn(ii+1) ! longitude boundaries of the physics grid
    71   REAL, INTENT (IN) :: airedyn(ii+1,jj+1) ! area of the dynamics grid (m2)
    72   REAL, INTENT (IN) :: cudyn((ii+1)*(jj+1)) ! cu coeff. (u_covariant = cu * u)
    73   REAL, INTENT (IN) :: cvdyn((ii+1)*jj) ! cv coeff. (v_covariant = cv * v)
    74   INTEGER, INTENT (IN) :: pdayref ! reference day of for the simulation
    75   REAL, INTENT (IN) :: ptimestep !physics time step (s)
    76   INTEGER, INTENT (IN) :: iflag_phys ! type of physics to be called
     56    REAL, INTENT (IN) :: prad ! radius of the planet (m)
     57    REAL, INTENT (IN) :: pg ! gravitational acceleration (m/s2)
     58    REAL, INTENT (IN) :: pr ! reduced gas constant R/mu
     59    REAL, INTENT (IN) :: pcpp ! specific heat Cp
     60    REAL, INTENT (IN) :: punjours ! length (in s) of a standard day
     61    INTEGER, INTENT (IN) :: nlayer ! number of atmospheric layers
     62    INTEGER, INTENT (IN) :: ii ! number of atmospheric columns along longitudes
     63    INTEGER, INTENT (IN) :: jj ! number of atompsheric columns along latitudes
     64    INTEGER, INTENT(IN) :: nbp ! number of physics columns for this MPI process
     65    INTEGER, INTENT(IN) :: communicator ! MPI communicator
     66    REAL, INTENT (IN) :: rlatudyn(jj + 1) ! latitudes of the physics grid
     67    REAL, INTENT (IN) :: rlatvdyn(jj) ! latitude boundaries of the physics grid
     68    REAL, INTENT (IN) :: rlonvdyn(ii + 1) ! longitudes of the physics grid
     69    REAL, INTENT (IN) :: rlonudyn(ii + 1) ! longitude boundaries of the physics grid
     70    REAL, INTENT (IN) :: airedyn(ii + 1, jj + 1) ! area of the dynamics grid (m2)
     71    REAL, INTENT (IN) :: cudyn((ii + 1) * (jj + 1)) ! cu coeff. (u_covariant = cu * u)
     72    REAL, INTENT (IN) :: cvdyn((ii + 1) * jj) ! cv coeff. (v_covariant = cv * v)
     73    INTEGER, INTENT (IN) :: pdayref ! reference day of for the simulation
     74    REAL, INTENT (IN) :: ptimestep !physics time step (s)
     75    INTEGER, INTENT (IN) :: iflag_phys ! type of physics to be called
    7776
    78   INTEGER :: ibegin, iend, offset
    79   INTEGER :: i,j,k
    80   CHARACTER (LEN=20) :: modname = 'iniphysiq'
    81   CHARACTER (LEN=80) :: abort_message
    82  
    83   LOGICAL :: slab_hdiff
    84   INTEGER :: slab_ekman
    85   CHARACTER (LEN = 6) :: type_ocean
     77    INTEGER :: ibegin, iend, offset
     78    INTEGER :: i, j, k
     79    CHARACTER (LEN = 20) :: modname = 'iniphysiq'
     80    CHARACTER (LEN = 80) :: abort_message
     81
     82    LOGICAL :: slab_hdiff
     83    INTEGER :: slab_ekman
     84    CHARACTER (LEN = 6) :: type_ocean
    8685
    8786#ifndef CPP_PARA
    88   INTEGER,PARAMETER :: mpi_rank=0
    89   INTEGER, PARAMETER :: mpi_size = 1
    90   INTEGER :: distrib_phys(mpi_rank:mpi_rank)=(jjm-1)*iim+2
     87    INTEGER, PARAMETER :: mpi_rank = 0
     88    INTEGER, PARAMETER :: mpi_size = 1
     89    INTEGER :: distrib_phys(mpi_rank:mpi_rank) = (jjm - 1) * iim + 2
    9190#endif
    9291
    93   ! --> initialize physics distribution, global fields and geometry
    94   ! (i.e. things in phy_common or dynphy_lonlat)
    95   CALL inigeomphy(ii,jj,nlayer, &
    96                nbp, communicator, &
    97                rlatudyn,rlatvdyn, &
    98                rlonudyn,rlonvdyn, &
    99                airedyn,cudyn,cvdyn)
     92    ! --> initialize physics distribution, global fields and geometry
     93    ! (i.e. things in phy_common or dynphy_lonlat)
     94    CALL inigeomphy(ii, jj, nlayer, &
     95            nbp, communicator, &
     96            rlatudyn, rlatvdyn, &
     97            rlonudyn, rlonvdyn, &
     98            airedyn, cudyn, cvdyn)
    10099
    101   ! --> now initialize things specific to the phylmd physics package
    102  
    103 !!$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/)
    104 !       Copy all threadprivate variables in temps_mod
    105 !$OMP PARALLEL DEFAULT(SHARED) COPYIN(annee_ref,day_ini,day_ref,start_time)
     100    ! --> now initialize things specific to the phylmd physics package
    106101
    107   ! Initialize physical constants in physics:
    108   CALL inifis(punjours,prad,pg,pr,pcpp)
     102    !!$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/)
     103    !   Copy all threadprivate variables in temps_mod
     104    !$OMP PARALLEL DEFAULT(SHARED) COPYIN(annee_ref,day_ini,day_ref,start_time)
    109105
    110   CALL init_time(annee_ref,day_ref,day_ini,start_time,nday,ptimestep)
     106    ! Initialize physical constants in physics:
     107    CALL inifis(punjours, prad, pg, pr, pcpp)
    111108
    112   ! Initialize dimphy module (unless in 1D where it has already been done)
    113 !  IF (klon_glo>1) CALL Init_dimphy(klon_omp,nlayer)
     109    CALL init_time(annee_ref, day_ref, day_ini, start_time, nday, ptimestep)
    114110
    115   ! Copy over "offline" settings
    116   CALL init_phystokenc(offline,istphy)
     111    ! Initialize dimphy module (unless in 1D where it has already been done)
     112    !  IF (klon_glo>1) CALL Init_dimphy(klon_omp,nlayer)
    117113
    118   ! Initialization for slab heat transport
    119   type_ocean="force"
    120   CALL getin_p('type_ocean',type_ocean)
    121   slab_hdiff=.FALSE.
    122   CALL getin_p('slab_hdiff',slab_hdiff)
    123   slab_ekman=0
    124   CALL getin_p('slab_ekman',slab_ekman)
    125   IF ((type_ocean=='slab').AND.(slab_hdiff.OR.(slab_ekman>0))) THEN
    126      CALL ini_slab_transp_geom(ip1jm,ip1jmp1,unsairez,fext,unsaire,&
    127                                   cu,cuvsurcv,cv,cvusurcu, &
    128                                   aire,apoln,apols, &
    129                                   aireu,airev,rlatvdyn,rad,omeg)
    130   END IF
     114    ! Copy over "offline" settings
     115    CALL init_phystokenc(offline, istphy)
    131116
    132   ! Initialize tracer names, numbers, etc. for physics
    133   CALL init_infotrac_phy
     117    ! Initialization for slab heat transport
     118    type_ocean = "force"
     119    CALL getin_p('type_ocean', type_ocean)
     120    slab_hdiff = .FALSE.
     121    CALL getin_p('slab_hdiff', slab_hdiff)
     122    slab_ekman = 0
     123    CALL getin_p('slab_ekman', slab_ekman)
     124    IF ((type_ocean=='slab').AND.(slab_hdiff.OR.(slab_ekman>0))) THEN
     125      CALL ini_slab_transp_geom(ip1jm, ip1jmp1, unsairez, fext, unsaire, &
     126              cu, cuvsurcv, cv, cvusurcu, &
     127              aire, apoln, apols, &
     128              aireu, airev, rlatvdyn, rad, omeg)
     129    END IF
    134130
    135   ! Initializations for Reprobus
    136   IF (type_trac == 'repr') THEN
     131    ! Initialize tracer names, numbers, etc. for physics
     132    CALL init_infotrac_phy
     133
     134    ! Initializations for Reprobus
     135    IF (type_trac == 'repr') THEN
    137136#ifdef REPROBUS
    138137    CALL Init_chem_rep_phys(klon_omp,nlayer)
    … …  
    141140          distrib_phys,communicator)
    142141#endif
    143   ENDIF
    144 !$OMP END PARALLEL
     142    ENDIF
     143    !$OMP END PARALLEL
    145144
    146 
    147   IF (type_trac == 'repr') THEN
     145    IF (type_trac == 'repr') THEN
    148146#ifdef REPROBUS
    149147    CALL init_reprobus_para( &
    … …  
    151149          distrib_phys,communicator)
    152150#endif
    153   ENDIF
     151    ENDIF
    154152
    155 !!$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/)
    156 !$OMP PARALLEL DEFAULT(SHARED)
    157   ! Additional initializations for aquaplanets
    158   IF (iflag_phys>=100) THEN
    159     CALL iniaqua(klon_omp,year_len,iflag_phys)
    160   END IF
     153    !!$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/)
     154    !$OMP PARALLEL DEFAULT(SHARED)
     155    ! Additional initializations for aquaplanets
     156    IF (iflag_phys>=100) THEN
     157      CALL iniaqua(klon_omp, year_len, iflag_phys)
     158    END IF
    161159
    162   IF (ANY(type_trac == ['inca','inco'])) THEN
    163      CALL init_inca_dim_reg(nbp_lon, nbp_lat - 1, &
    164           rlonudyn, rlatudyn, rlonvdyn, rlatvdyn)
    165   END IF
     160    IF (ANY(type_trac == ['inca', 'inco'])) THEN
     161      CALL init_inca_dim_reg(nbp_lon, nbp_lat - 1, &
     162              rlonudyn, rlatudyn, rlonvdyn, rlatvdyn)
     163    END IF
    166164
    167 !$OMP END PARALLEL
     165    !$OMP END PARALLEL
    168166
    169 END SUBROUTINE iniphysiq
     167  END SUBROUTINE iniphysiq
    170168
    171169END MODULE iniphysiq_mod

  • LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/phylmd/init_ssrf_m.F90

    r5112 r5118  
    11MODULE init_ssrf_m
    22
    3 !*******************************************************************************
     3  !*******************************************************************************
    44
    5   USE indice_sol_mod,     ONLY: is_ter, is_oce, is_oce, is_lic, epsfra
    6   USE dimphy,             ONLY: klon, zmasq
     5  USE indice_sol_mod, ONLY: is_ter, is_oce, is_oce, is_lic, epsfra
     6  USE dimphy, ONLY: klon, zmasq
    77  USE phys_state_var_mod, ONLY: pctsrf
    8   USE lmdz_geometry,       ONLY: longitude_deg, latitude_deg
    9   USE grid_atob_m,        ONLY: grille_m
    10   USE ioipsl,             ONLY: flininfo, flinopen, flinget, flinclo
     8  USE lmdz_geometry, ONLY: longitude_deg, latitude_deg
     9  USE grid_atob_m, ONLY: grille_m
     10  USE ioipsl, ONLY: flininfo, flinopen, flinget, flinclo
    1111  USE lmdz_ioipsl_getin_p, ONLY: getin_p
    12   USE comconst_mod,       ONLY: im, pi
    13   USE surface_data,       ONLY: landice_opt
     12  USE comconst_mod, ONLY: im, pi
     13  USE surface_data, ONLY: landice_opt
     14  USE lmdz_iniprint, ONLY: lunout, prt_level
    1415
    15   CHARACTER(LEN=256), PARAMETER :: icefname="landiceref.nc", icevar="landice"
     16  CHARACTER(LEN = 256), PARAMETER :: icefname = "landiceref.nc", icevar = "landice"
    1617  PRIVATE
    1718  PUBLIC :: start_init_subsurf
    18   include "iniprint.h"
    1919  include "dimensions.h"
    2020  include "paramet.h"
    … …  
    2323CONTAINS
    2424
    25 !-------------------------------------------------------------------------------
     25  !-------------------------------------------------------------------------------
    2626
    27 SUBROUTINE start_init_subsurf(known_mask)
     27  SUBROUTINE start_init_subsurf(known_mask)
    2828
    29 !-------------------------------------------------------------------------------
    30 ! Purpose: Subsurfaces initialization.
    31 !-------------------------------------------------------------------------------
    32 ! Comment: Called by etat0phys_netcdf ; also called by limit_netcdf in case
    33 !          no starting states are required (ok_etat0==.FALSE.).
    34 !-------------------------------------------------------------------------------
    35   IMPLICIT NONE
    36 !-------------------------------------------------------------------------------
    37 ! Arguments:
    38   LOGICAL, INTENT(IN) :: known_mask
    39 !-------------------------------------------------------------------------------
    40 ! Local variables:
    41   INTEGER          :: iml_lic, jml_lic
    42   INTEGER          :: fid, llm_tmp, ttm_tmp, itaul(1), ji, j
    43   REAL, ALLOCATABLE :: dlon_lic(:), lon_lic(:,:), fraclic (:,:)
    44   REAL, ALLOCATABLE :: dlat_lic(:), lat_lic(:,:), flic_tmp(:,:), vtmp(:,:)
    45   REAL              :: date, lev(1), dt, deg2rad
    46   LOGICAL          :: no_ter_antartique   ! If true, no land points are allowed at Antartic
    47 !-------------------------------------------------------------------------------
    48   deg2rad= pi/180.0
     29    !-------------------------------------------------------------------------------
     30    ! Purpose: Subsurfaces initialization.
     31    !-------------------------------------------------------------------------------
     32    ! Comment: Called by etat0phys_netcdf ; also called by limit_netcdf in case
     33    !          no starting states are required (ok_etat0==.FALSE.).
     34    !-------------------------------------------------------------------------------
     35    IMPLICIT NONE
     36    !-------------------------------------------------------------------------------
     37    ! Arguments:
     38    LOGICAL, INTENT(IN) :: known_mask
     39    !-------------------------------------------------------------------------------
     40    ! Local variables:
     41    INTEGER :: iml_lic, jml_lic
     42    INTEGER :: fid, llm_tmp, ttm_tmp, itaul(1), ji, j
     43    REAL, ALLOCATABLE :: dlon_lic(:), lon_lic(:, :), fraclic (:, :)
     44    REAL, ALLOCATABLE :: dlat_lic(:), lat_lic(:, :), flic_tmp(:, :), vtmp(:, :)
     45    REAL :: date, lev(1), dt, deg2rad
     46    LOGICAL :: no_ter_antartique   ! If true, no land points are allowed at Antartic
     47    !-------------------------------------------------------------------------------
     48    deg2rad = pi / 180.0
    4949
    50 !--- Physical grid points coordinates
    51   DO j=2,jjm; latitude_deg((j-2)*iim+2:(j-1)*iim+1)=rlatu(j);    END DO
    52   DO j=2,jjm; longitude_deg((j-2)*iim+2:(j-1)*iim+1)=rlonv(1:im); END DO
    53   latitude_deg(1) = pi/2.; latitude_deg(klon) = - pi/2.
    54   latitude_deg(:)=latitude_deg(:)/deg2rad
    55   longitude_deg(1) = 0.0;   longitude_deg(klon) = 0.0;
    56   longitude_deg(:)=longitude_deg(:)/deg2rad
     50    !--- Physical grid points coordinates
     51    DO j = 2, jjm; latitude_deg((j - 2) * iim + 2:(j - 1) * iim + 1) = rlatu(j);
     52    END DO
     53    DO j = 2, jjm; longitude_deg((j - 2) * iim + 2:(j - 1) * iim + 1) = rlonv(1:im);
     54    END DO
     55    latitude_deg(1) = pi / 2.; latitude_deg(klon) = - pi / 2.
     56    latitude_deg(:) = latitude_deg(:) / deg2rad
     57    longitude_deg(1) = 0.0;   longitude_deg(klon) = 0.0;
     58    longitude_deg(:) = longitude_deg(:) / deg2rad
    5759
    58 ! Compute ground geopotential, sub-cells quantities and possibly the mask.
    59 ! Sub-surfaces initialization
    60 !*******************************************************************************
    61   IF (landice_opt < 2) THEN
    62      ! Continue with reading landice.nc file
    63      WRITE(lunout,*)"Read landice.nc file to attribute is_lic fraction"
     60    ! Compute ground geopotential, sub-cells quantities and possibly the mask.
     61    ! Sub-surfaces initialization
     62    !*******************************************************************************
     63    IF (landice_opt < 2) THEN
     64      ! Continue with reading landice.nc file
     65      WRITE(lunout, *)"Read landice.nc file to attribute is_lic fraction"
    6466
    65      !--- Read and interpolate on model T-grid soil fraction and soil ice fraction.
    66      CALL flininfo(icefname, iml_lic, jml_lic, llm_tmp, ttm_tmp, fid)
    67      ALLOCATE(lat_lic(iml_lic,jml_lic),lon_lic(iml_lic,jml_lic))
    68      ALLOCATE(fraclic(iml_lic,jml_lic))
    69      CALL flinopen(icefname, .FALSE., iml_lic, jml_lic, llm_tmp, &
    70                           lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
    71      CALL flinget(fid, icevar, iml_lic, jml_lic, llm_tmp, ttm_tmp, 1,1, fraclic)
    72      CALL flinclo(fid)
    73      WRITE(lunout,*)'landice dimensions: iml_lic, jml_lic : ',iml_lic,jml_lic
    74      
    75      ALLOCATE(dlon_lic(iml_lic),dlat_lic(jml_lic))
    76      dlon_lic(:)=lon_lic(:,1); IF(MAXVAL(dlon_lic)>pi) dlon_lic=dlon_lic*pi/180.
    77      dlat_lic(:)=lat_lic(1,:); IF(MAXVAL(dlat_lic)>pi) dlat_lic=dlat_lic*pi/180.
    78      DEALLOCATE(lon_lic,lat_lic); ALLOCATE(flic_tmp(iip1,jjp1))
    79      CALL grille_m(dlon_lic,dlat_lic,fraclic,rlonv(1:iim),rlatu,flic_tmp(1:iim,:))
    80      flic_tmp(iip1,:)=flic_tmp(1,:)
    81      
    82      !--- To the physical grid
    83      pctsrf(:,:) = 0.
    84      CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp, pctsrf(:,is_lic))
    85      DEALLOCATE(flic_tmp)
    86      
    87      !--- Adequation with soil/sea mask
    88      WHERE(pctsrf(:,is_lic)<EPSFRA) pctsrf(:,is_lic)=0.
    89      WHERE(zmasq(:)<EPSFRA)         pctsrf(:,is_lic)=0.
    90      pctsrf(:,is_ter)=zmasq(:)
    91      DO ji=1,klon
    92         IF(zmasq(ji)>EPSFRA) THEN 
    93            IF(pctsrf(ji,is_lic)>=zmasq(ji)) THEN
    94               pctsrf(ji,is_lic)=zmasq(ji)
    95               pctsrf(ji,is_ter)=0.
    96            ELSE
    97               pctsrf(ji,is_ter)=zmasq(ji)-pctsrf(ji,is_lic)
    98               IF(pctsrf(ji,is_ter)<EPSFRA) THEN
    99                  pctsrf(ji,is_ter)=0.
    100                  pctsrf(ji,is_lic)=zmasq(ji)
    101               END IF
    102            END IF
     67      !--- Read and interpolate on model T-grid soil fraction and soil ice fraction.
     68      CALL flininfo(icefname, iml_lic, jml_lic, llm_tmp, ttm_tmp, fid)
     69      ALLOCATE(lat_lic(iml_lic, jml_lic), lon_lic(iml_lic, jml_lic))
     70      ALLOCATE(fraclic(iml_lic, jml_lic))
     71      CALL flinopen(icefname, .FALSE., iml_lic, jml_lic, llm_tmp, &
     72              lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
     73      CALL flinget(fid, icevar, iml_lic, jml_lic, llm_tmp, ttm_tmp, 1, 1, fraclic)
     74      CALL flinclo(fid)
     75      WRITE(lunout, *)'landice dimensions: iml_lic, jml_lic : ', iml_lic, jml_lic
     76
     77      ALLOCATE(dlon_lic(iml_lic), dlat_lic(jml_lic))
     78      dlon_lic(:) = lon_lic(:, 1); IF(MAXVAL(dlon_lic)>pi) dlon_lic = dlon_lic * pi / 180.
     79      dlat_lic(:) = lat_lic(1, :); IF(MAXVAL(dlat_lic)>pi) dlat_lic = dlat_lic * pi / 180.
     80      DEALLOCATE(lon_lic, lat_lic); ALLOCATE(flic_tmp(iip1, jjp1))
     81      CALL grille_m(dlon_lic, dlat_lic, fraclic, rlonv(1:iim), rlatu, flic_tmp(1:iim, :))
     82      flic_tmp(iip1, :) = flic_tmp(1, :)
     83
     84      !--- To the physical grid
     85      pctsrf(:, :) = 0.
     86      CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp, pctsrf(:, is_lic))
     87      DEALLOCATE(flic_tmp)
     88
     89      !--- Adequation with soil/sea mask
     90      WHERE(pctsrf(:, is_lic)<EPSFRA) pctsrf(:, is_lic) = 0.
     91      WHERE(zmasq(:)<EPSFRA)         pctsrf(:, is_lic) = 0.
     92      pctsrf(:, is_ter) = zmasq(:)
     93      DO ji = 1, klon
     94        IF(zmasq(ji)>EPSFRA) THEN
     95          IF(pctsrf(ji, is_lic)>=zmasq(ji)) THEN
     96            pctsrf(ji, is_lic) = zmasq(ji)
     97            pctsrf(ji, is_ter) = 0.
     98          ELSE
     99            pctsrf(ji, is_ter) = zmasq(ji) - pctsrf(ji, is_lic)
     100            IF(pctsrf(ji, is_ter)<EPSFRA) THEN
     101              pctsrf(ji, is_ter) = 0.
     102              pctsrf(ji, is_lic) = zmasq(ji)
     103            END IF
     104          END IF
    103105        END IF
    104      END DO
    105      
    106      
    107      !--- Option no_ter_antartique removes all land fractions souther than 60S.
    108      !--- Land ice is set instead of the land fractions on these latitudes.
    109      !--- The ocean and sea-ice fractions are not changed.
    110      no_ter_antartique=.FALSE.
    111      CALL getin_p('no_ter_antartique',no_ter_antartique)
    112      WRITE(lunout,*)"no_ter_antartique=",no_ter_antartique
    113      IF (no_ter_antartique) THEN
     106      END DO
     107
     108
     109      !--- Option no_ter_antartique removes all land fractions souther than 60S.
     110      !--- Land ice is set instead of the land fractions on these latitudes.
     111      !--- The ocean and sea-ice fractions are not changed.
     112      no_ter_antartique = .FALSE.
     113      CALL getin_p('no_ter_antartique', no_ter_antartique)
     114      WRITE(lunout, *)"no_ter_antartique=", no_ter_antartique
     115      IF (no_ter_antartique) THEN
    114116        ! Remove all land fractions souther than 60S and set land-ice instead
    115         WRITE(lunout,*) "Remove land fractions souther than 60deg south by increasing"
    116         WRITE(lunout,*) "the continental ice fractions. No land can now be found at Antartic."
    117         DO ji=1, klon
    118            IF (latitude_deg(ji)<-60.0) THEN
    119               pctsrf(ji,is_lic) = pctsrf(ji,is_lic) + pctsrf(ji,is_ter)
    120               pctsrf(ji,is_ter) = 0
    121            END IF
     117        WRITE(lunout, *) "Remove land fractions souther than 60deg south by increasing"
     118        WRITE(lunout, *) "the continental ice fractions. No land can now be found at Antartic."
     119        DO ji = 1, klon
     120          IF (latitude_deg(ji)<-60.0) THEN
     121            pctsrf(ji, is_lic) = pctsrf(ji, is_lic) + pctsrf(ji, is_ter)
     122            pctsrf(ji, is_ter) = 0
     123          END IF
    122124        END DO
    123      END IF
    124      
    125   ELSE
    126      ! landice_opt=2 and higher
    127      WRITE(lunout,*) 'No landice is attributed is_lic sub-surface because landice_opt=2 or higher.'
    128      WRITE(lunout,*) 'This means that the land model will handel land ice as well as all other land areas.'
    129      pctsrf(:,is_ter) = zmasq(:)
    130      pctsrf(:,is_lic) = 0.0
    131   END IF
     125      END IF
    132126
    133 !--- Sub-surface ocean and sea ice (sea ice set to zero for start).
    134   pctsrf(:,is_oce)=(1.-zmasq(:))
    135   WHERE(pctsrf(:,is_oce)<EPSFRA) pctsrf(:,is_oce)=0.
    136   IF(known_mask) pctsrf(:,is_oce)=1-zmasq(:)
     127    ELSE
     128      ! landice_opt=2 and higher
     129      WRITE(lunout, *) 'No landice is attributed is_lic sub-surface because landice_opt=2 or higher.'
     130      WRITE(lunout, *) 'This means that the land model will handel land ice as well as all other land areas.'
     131      pctsrf(:, is_ter) = zmasq(:)
     132      pctsrf(:, is_lic) = 0.0
     133    END IF
    137134
    138 !--- It is checked that the sub-surfaces sum is equal to 1.
    139   ji=COUNT((ABS(SUM(pctsrf(:,:),dim=2))-1.0)>EPSFRA)
    140   IF(ji/=0) WRITE(lunout,*) 'Sub-cell distribution problem for ',ji,' points'
     135    !--- Sub-surface ocean and sea ice (sea ice set to zero for start).
     136    pctsrf(:, is_oce) = (1. - zmasq(:))
     137    WHERE(pctsrf(:, is_oce)<EPSFRA) pctsrf(:, is_oce) = 0.
     138    IF(known_mask) pctsrf(:, is_oce) = 1 - zmasq(:)
    141139
    142 END SUBROUTINE start_init_subsurf
     140    !--- It is checked that the sub-surfaces sum is equal to 1.
     141    ji = COUNT((ABS(SUM(pctsrf(:, :), dim = 2)) - 1.0)>EPSFRA)
     142    IF(ji/=0) WRITE(lunout, *) 'Sub-cell distribution problem for ', ji, ' points'
    143143
    144 !-------------------------------------------------------------------------------
     144  END SUBROUTINE start_init_subsurf
     145
     146  !-------------------------------------------------------------------------------
    145147
    146148END MODULE init_ssrf_m

  • LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/phylmd/limit_netcdf.f90

    r5117 r5118  
    8181    USE phys_cal_mod, ONLY: calend
    8282    USE lmdz_abort_physic, ONLY: abort_physic
     83    USE lmdz_iniprint, ONLY: lunout, prt_level
    8384    IMPLICIT NONE
    8485    !-------------------------------------------------------------------------------
    8586    ! Arguments:
    86     include "iniprint.h"
    8787    include "dimensions.h"
    8888    include "paramet.h"
Note: See TracChangeset for help on using the changeset viewer.