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

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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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
Note: See TracChangeset for help on using the changeset viewer.