Changeset 11586


Ignore:
Timestamp:
2019-09-20T17:28:02+02:00 (14 months ago)
Author:
gsamson
Message:

dev_r11265_ABL : see #2131

  • merge HPC-13_IRRMANN_BDY_optimization branch @ r11535 (last commit) with dev_r11265_ABL branch @ r11414 (except doc directory)
  • change ORCA2 results due to changes in HPC-13_IRRMANN_BDY_optimization branch
Location:
NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D
Files:
2 deleted
89 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/ORCA2_ICE_PISCES/EXPREF/context_nemo.xml

    r11413 r11586  
    4242      <axis id="iax_20C" long_name="20 degC isotherm"  unit="degC"              /> 
    4343      <axis id="iax_28C" long_name="28 degC isotherm"  unit="degC"              /> 
    44       <!-- ABL vertical axis definition --> 
    45       <axis id="ght_abl" long_name="ABL Vertical T levels" unit="m" positive="up"   /> 
    46       <axis id="ghw_abl" long_name="ABL Vertical W levels" unit="m" positive="up"   /> 
    4744    </axis_definition> 
    4845  
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/field_def_nemo-ice.xml

    r11413 r11586  
    4646      
    4747     <!-- melt ponds --> 
    48      <field id="iceapnd"      long_name="melt pond fraction"                                      standard_name="sea_ice_meltpond_fraction"                 unit="%" />  
     48     <field id="iceapnd"      long_name="melt pond concentration"                                 standard_name="sea_ice_meltpond_concentration"            unit=""  />  
     49          <field id="icehpnd"      long_name="melt pond depth"                                         standard_name="sea_ice_meltpond_depth"                    unit="m" />  
    4950          <field id="icevpnd"      long_name="melt pond volume"                                        standard_name="sea_ice_meltpond_volume"                   unit="m" />  
    5051      
     
    328329     <field field_ref="icesalt"          name="sisali" /> 
    329330     <field field_ref="iceapnd"          name="siapnd" /> 
     331     <field field_ref="icehpnd"          name="sihpnd" /> 
    330332     <field field_ref="icevpnd"          name="sivpnd" /> 
    331333          <field field_ref="iceage"           name="siage"  /> 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/field_def_nemo-oce.xml

    r11413 r11586  
    237237       
    238238      <!-- SBC --> 
     239 
    239240      <field_group id="SBC" > <!-- time step automaticaly defined based on nn_fsbc --> 
    240241 
     
    300301          <field id="iceshelf_cea" long_name="Iceshelf"                                  standard_name="water_flux_into_sea_water_from_iceshelf"   unit="kg/m2/s"  /> 
    301302 
    302  
    303303          <!-- available if key_oasis3 + conservative method --> 
    304304          <field id="rain"          long_name="Liquid precipitation"                                     standard_name="rainfall_flux"                                                                 unit="kg/m2/s"  /> 
     
    348348      <!-- ABL --> 
    349349      <field_group id="ABL" > <!-- time step automaticaly defined based on nn_fsbc --> 
     350 
    350351   <!-- variables available with ABL on atmospheric T grid--> 
    351352   <field_group id="grid_ABL3D" grid_ref="grid_TA_3D" > 
     
    354355          <field id="t_abl"      long_name="ABL potential temperature"     standard_name="abl_theta"      unit="K"        /> 
    355356          <field id="q_abl"      long_name="ABL specific humidity"         standard_name="abl_qspe"       unit="kg/kg"    /> 
     357          <!-- debug (to be removed) --> 
     358          <field id="u_dta"      long_name="DTA i-horizontal velocity"     standard_name="dta_x_velocity" unit="m/s"      /> 
     359          <field id="v_dta"      long_name="DTA j-horizontal velocity"     standard_name="dta_y_velocity" unit="m/s"      /> 
     360          <field id="t_dta"      long_name="DTA potential temperature"     standard_name="dta_theta"      unit="K"        /> 
     361          <field id="q_dta"      long_name="DTA specific humidity"         standard_name="dta_qspe"       unit="kg/kg"    /> 
    356362          <field id="coeft"      long_name="ABL nudging coefficient"       standard_name="coeft"          unit=""         /> 
    357363   </field_group> 
     364 
    358365   <field_group id="grid_ABL2D" grid_ref="grid_TA_2D" > 
    359366          <field id="pblh"       long_name="ABL height"                    standard_name="abl_height"     unit="m"        /> 
     367          <field id="uz1_abl"    long_name="ABL i-horizontal velocity"     standard_name="abl_x_velocity" unit="m/s"      /> 
     368          <field id="vz1_abl"    long_name="ABL j-horizontal velocity"     standard_name="abl_y_velocity" unit="m/s"      /> 
     369          <field id="uvz1_abl"   long_name="ABL wind speed module"         standard_name="abl_wind_speed" unit="m/s"       > sqrt( uz1_abl^2 + vz1_abl^2 ) </field> 
     370          <field id="tz1_abl"    long_name="ABL potential temperature"     standard_name="abl_theta"      unit="K"        /> 
     371          <field id="qz1_abl"    long_name="ABL specific humidity"         standard_name="abl_qspe"       unit="kg/kg"    /> 
     372          <field id="uz1_dta"    long_name="DTA i-horizontal velocity"     standard_name="dta_x_velocity" unit="m/s"      /> 
     373          <field id="vz1_dta"    long_name="DTA j-horizontal velocity"     standard_name="dta_y_velocity" unit="m/s"      /> 
     374          <field id="uvz1_dta"   long_name="DTA wind speed module"         standard_name="dta_wind_speed" unit="m/s"       > sqrt( uz1_dta^2 + vz1_dta^2 ) </field>  
     375          <field id="tz1_dta"    long_name="DTA potential temperature"     standard_name="dta_theta"      unit="K"        /> 
     376          <field id="qz1_dta"    long_name="DTA specific humidity"         standard_name="dta_qspe"       unit="kg/kg"    /> 
     377          <!-- debug (to be removed) --> 
     378          <field id="uz1_geo"    long_name="GEO i-horizontal velocity"     standard_name="geo_x_velocity" unit="m/s"      /> 
     379          <field id="vz1_geo"    long_name="GEO j-horizontal velocity"     standard_name="geo_y_velocity" unit="m/s"      /> 
     380          <field id="uvz1_geo"   long_name="GEO wind speed module"         standard_name="geo_wind_speed" unit="m/s"       > sqrt( uz1_geo^2 + vz1_geo^2 ) </field>  
    360381   </field_group> 
     382 
    361383      </field_group> <!-- ABL --> 
     384 
    362385       
    363386      <!-- U grid --> 
     
    375398        <field id="uocet"        long_name="ocean transport along i-axis times temperature (CRS)"                                               unit="degC*m/s"   grid_ref="grid_U_3D" /> 
    376399        <field id="uoces"        long_name="ocean transport along i-axis times salinity (CRS)"                                                  unit="1e-3*m/s"   grid_ref="grid_U_3D" /> 
     400        <field id="ssuww"        long_name="ocean surface wind work along i-axis"                   standard_name="surface_x_wind_work"         unit="N/m*s"                            > utau * ssu </field> 
    377401 
    378402        <!-- u-eddy diffusivity coefficients (available if ln_traldf_OFF=F) --> 
     
    429453        <field id="vocet"        long_name="ocean transport along j-axis times temperature (CRS)"                                               unit="degC*m/s"   grid_ref="grid_V_3D" /> 
    430454        <field id="voces"        long_name="ocean transport along j-axis times salinity (CRS)"                                                  unit="1e-3*m/s"   grid_ref="grid_V_3D" /> 
     455        <field id="ssvww"        long_name="ocean surface wind work along j-axis"                   standard_name="surface_y_wind_work"         unit="N/m*s"                            > vtau * ssv </field> 
     456 
    431457 
    432458        <!-- v-eddy diffusivity coefficients (available if ln_traldf_OFF=F) --> 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/namelist_ice_ref

    r11413 r11586  
    7070   ln_str_H79       = .true.          !  ice strength param.: Hibler_79   => P = pstar*<h>*exp(-c_rhg*A)                       
    7171      rn_pstar      =   2.0e+04       !     ice strength thickness parameter [N/m2] 
    72       rn_crhg       =   20.0          !     ice strength conc. parameter (-) 
     72      rn_crhg       =  20.0           !     ice strength conc. parameter (-) 
    7373                   ! -- ice_rdgrft -- ! 
    7474   rn_csrdg         =   0.5           !  fraction of shearing energy contributing to ridging 
     
    177177&namthd_pnd     !   Melt ponds 
    178178!------------------------------------------------------------------------------ 
    179    ln_pnd_H12       = .false.         !  activate evolutive melt ponds (from Holland et al 2012) 
    180    ln_pnd_CST       = .false.         !  activate constant  melt ponds 
    181       rn_apnd       =   0.2           !     prescribed pond fraction, at Tsu=0 degC 
    182       rn_hpnd       =   0.05          !     prescribed pond depth,    at Tsu=0 degC 
    183    ln_pnd_alb       = .false.         !  melt ponds affect albedo or not 
     179   ln_pnd           = .false.         !  activate melt ponds or not 
     180     ln_pnd_H12     = .false.         !  activate evolutive melt ponds (from Holland et al 2012) 
     181     ln_pnd_CST     = .false.         !  activate constant  melt ponds 
     182       rn_apnd      =   0.2           !     prescribed pond fraction, at Tsu=0 degC 
     183       rn_hpnd      =   0.05          !     prescribed pond depth,    at Tsu=0 degC 
     184     ln_pnd_alb     = .false.         !  melt ponds affect albedo or not 
    184185/ 
    185186!------------------------------------------------------------------------------ 
     
    207208   rn_hpd_ini_n     =   0.05          !  initial pond depth          (m), North 
    208209   rn_hpd_ini_s     =   0.05          !        "            "             South 
    209    !        ! if ln_iceini_file=T 
     210   ! -- for ln_iceini_file = T 
    210211   sn_hti = 'Ice_initialization'    , -12 ,'hti'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    211212   sn_hts = 'Ice_initialization'    , -12 ,'hts'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
     
    215216   sn_tsu = 'Ice_initialization'    , -12 ,'tsu'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    216217   sn_tms = 'NOT USED'              , -12 ,'tms'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
     218   !      melt ponds (be careful, sn_apd is the pond concentration (not fraction), so it differs from rn_apd) 
    217219   sn_apd = 'NOT USED'              , -12 ,'apd'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    218220   sn_hpd = 'NOT USED'              , -12 ,'hpd'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
     
    232234&namdia         !   Diagnostics 
    233235!------------------------------------------------------------------------------ 
    234    ln_icediachk     = .false.         !  check online the heat, mass & salt budgets (T) or not (F) 
     236   ln_icediachk     = .false.         !  check online the heat, mass & salt budgets at each time step 
     237      !                               !     rate of ice spuriously gained/lost. For ex., rn_icechk=1. <=> 1mm/year, rn_icechk=0.1 <=> 1mm/10years                                    
     238      rn_icechk_cel =  1.             !     check at any gridcell           => stops the code if violated (and writes a file) 
     239      rn_icechk_glo =  0.1            !     check over the entire ice cover => only prints warnings 
    235240   ln_icediahsb     = .false.         !  output the heat, mass & salt budgets (T) or not (F) 
    236241   ln_icectl        = .false.         !  ice points output for debug (T or F) 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SHARED/namelist_ref

    r11413 r11586  
    199199   ln_flx      = .false.   !  flux formulation                          (T => fill namsbc_flx ) 
    200200   ln_blk      = .false.   !  Bulk formulation                          (T => fill namsbc_blk ) 
     201   ln_abl      = .false.   !  ABL  formulation                          (T => fill namsbc_abl ) 
    201202      !              ! Type of coupling (Ocean/Ice/Atmosphere) : 
    202203   ln_cpl      = .false.   !  atmosphere coupled   formulation          ( requires key_oasis3 ) 
     
    280281   sn_snow     = 'ncar_precip.15JUNE2009_fill',   -1.        , 'SNOW'    ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , '' 
    281282   sn_slp      = 'slp.15JUNE2009_fill'        ,    6.        , 'SLP'     ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , '' 
     283/ 
     284!----------------------------------------------------------------------- 
     285&namsbc_abl    !   Atmospheric Boundary Layer formulation           (ln_abl = T) 
     286!----------------------------------------------------------------------- 
     287   cn_dir         = './'      !  root directory for the location of the ABL grid file 
     288   cn_dom         = 'dom_cfg_abl.nc' 
     289   ln_hpgls_frc   = .false. 
     290   ln_geos_winds  = .false. 
     291   nn_dyn_restore = 2         ! restoring option for dynamical ABL variables: = 0 no restoring 
     292                              !                                               = 1 equatorial restoring 
     293                              !                                               = 2 global restoring 
     294   rn_ldyn_min   =  4.5       !  magnitude of the nudging on ABL dynamics at the bottom of the ABL   [hour] 
     295   rn_ldyn_max   =  1.5       !  magnitude of the nudging on ABL dynamics at the top of the ABL   [hour] 
     296   rn_ltra_min   =  4.5       !  magnitude of the nudging on ABL tracers  at the bottom of the ABL   [hour] 
     297   rn_ltra_max   =  1.5       !  magnitude of the nudging on ABL tracers  at the top of the ABL   [hour] 
     298   nn_amxl        = 0         ! mixing length: = 0 Deardorff 80 length-scale 
     299                              !                = 1 length-scale based on the distance to the PBL height 
     300                              !                = 2 Bougeault & Lacarrere 89 length-scale 
     301   rn_Cm         = 0.0667     ! 0.126 in MesoNH 
     302   rn_Ct         = 0.1667     ! 0.143 in MesoNH 
     303   rn_Ce         = 0.4        ! 0.4   in MesoNH 
     304   rn_Ceps       = 0.7        ! 0.85  in MesoNH 
     305   rn_Rod        = 0.15       ! c0 in RMCA17 mixing length formulation (not yet implemented) 
     306   rn_Ric        = 0.139      !  Critical Richardson number (to compute PBL height and diffusivities) 
    282307/ 
    283308!----------------------------------------------------------------------- 
     
    600625   nn_ice_dta    =  0         !  = 0, bdy data are equal to the initial state 
    601626   !                          !  = 1, bdy data are read in 'bdydata   .nc' files 
    602    rn_ice_tem    = 270.       !  si3 only: arbitrary temperature of incoming sea ice 
    603    rn_ice_sal    = 10.        !  si3 only:      --   salinity           -- 
    604    rn_ice_age    = 30.        !  si3 only:      --   age                -- 
    605627   ! 
    606628   ln_tra_dmp    =.false.     !  open boudaries conditions for tracers 
     
    632654   bn_sal      = 'amm12_bdyT_tra'        ,         24.       , 'vosaline',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
    633655!* for si3 
    634 !   bn_a_i     = 'amm12_bdyT_ice'        ,         24.       , 'ileadfra',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
    635 !   bn_h_i     = 'amm12_bdyT_ice'        ,         24.       , 'iicethic',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
    636 !   bn_h_s     = 'amm12_bdyT_ice'        ,         24.       , 'isnowthi',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
     656   bn_a_i      = 'amm12_bdyT_ice'        ,         24.       , 'siconc'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
     657   bn_h_i      = 'amm12_bdyT_ice'        ,         24.       , 'sithic'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
     658   bn_h_s      = 'amm12_bdyT_ice'        ,         24.       , 'snthic'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
     659   bn_t_i      = 'NOT USED'              ,         24.       , 'sitemp'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
     660   bn_t_s      = 'NOT USED'              ,         24.       , 'sntemp'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
     661   bn_tsu      = 'NOT USED'              ,         24.       , 'sittop'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
     662   bn_s_i      = 'NOT USED'              ,         24.       , 'sisalt'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
     663   ! melt ponds (be careful, bn_aip is the pond concentration (not fraction), so it differs from rn_iceapnd) 
     664   bn_aip      = 'NOT USED'              ,         24.       , 'siapnd'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
     665   bn_hip      = 'NOT USED'              ,         24.       , 'sihpnd'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
     666   ! if bn_t_i etc are "not used", then define arbitrary temperatures and salinity and ponds 
     667   rn_ice_tem  = 270.         !  arbitrary temperature               of incoming sea ice 
     668   rn_ice_sal  = 10.          !       --   salinity                            -- 
     669   rn_ice_age  = 30.          !       --   age                                 -- 
     670   rn_ice_apnd = 0.2          !       --   pond fraction = a_ip/a_i            -- 
     671   rn_ice_hpnd = 0.05         !       --   pond depth                          -- 
    637672/ 
    638673!----------------------------------------------------------------------- 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SPITZ12/EXPREF/namelist_cfg

    r11413 r11586  
    170170   nn_ice_dta    =  1         !  = 0, bdy data are equal to the initial state 
    171171   !                          !  = 1, bdy data are read in 'bdydata   .nc' files 
    172    rn_ice_tem    = 267.       !  si3 only: arbitrary temperature of incoming sea ice 
    173    rn_ice_sal    =   6.       !  si3 only:      --   salinity           -- 
    174    rn_ice_age    = 365.       !  si3 only:      --   age                -- 
    175    ! 
    176172   nn_rimwidth   = 1          !  width of the relaxation zone 
    177173   ln_vol        = .false.    !  total volume correction (see nn_volctl parameter) 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/SPITZ12/EXPREF/namelist_ice_cfg

    r11263 r11586  
    8080&namthd_pnd     !   Melt ponds 
    8181!------------------------------------------------------------------------------ 
    82    ln_pnd_H12       = .true.          !  activate evolutive melt ponds (from Holland et al 2012) 
    83    ln_pnd_alb       = .true.          !  melt ponds affect albedo or not 
     82   ln_pnd           = .true.          !  activate melt ponds or not 
     83     ln_pnd_H12     = .true.          !  activate evolutive melt ponds (from Holland et al 2012) 
     84     ln_pnd_alb     = .true.          !  melt ponds affect albedo or not 
    8485/ 
    8586 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/cfgs/ref_cfgs.txt

    r11413 r11586  
    1414eORCA025_ICE OCE ICE 
    1515eORCA025_ICE_ABL OCE ICE ABL 
     16eORCA025_SAS_ICE_ABL OCE SAS ICE ABL 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/ablmod.F90

    r11363 r11586  
    563563            DO ji = 2, jpim1   
    564564                
    565             zztmp1 = 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) ) 
    566             zztmp2 = 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) 
     565               zztmp1 = 0.5_wp * ( u_abl(ji+1,jj,2,nt_a) + u_abl(ji,jj,2,nt_a) ) 
     566               zztmp2 = 0.5_wp * ( v_abl(ji,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) ) 
    567567             
    568             ptaui_ice(ji,jj) = 0.5_wp * (  zrhoa(ji+1,jj) * pCd_du_ice(ji+1,jj)             & 
     568               ptaui_ice(ji,jj) = 0.5_wp * (  zrhoa(ji+1,jj) * pCd_du_ice(ji+1,jj)             & 
    569569                  &                      +    zrhoa(ji  ,jj) * pCd_du_ice(ji  ,jj)  )          & 
    570570                  &         * ( zztmp1 - rn_vfac * pssu_ice(ji,jj) ) 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ABL/sbcabl.F90

    r11363 r11586  
    214214            ! Check that restoring coefficients are between 0 and 1 
    215215            !IF( zcff1 > 1._wp .OR. zcff1 < 0._wp )   & 
    216             IF( zcff1 > nn_fsbc .OR. zcff1 < 0._wp )   & 
     216            !IF( zcff1 > nn_fsbc .OR. zcff1 < 0._wp )   & 
     217            IF( zcff1 - nn_fsbc > 0.001_wp .OR. zcff1 < 0._wp )   & 
    217218               &                   CALL ctl_stop( 'abl_init : wrong value for rn_ldyn_max' ) 
    218219            !IF( zcff  > 1._wp .OR. zcff  < 0._wp )   & 
    219             IF( zcff  > nn_fsbc .OR. zcff  < 0._wp )   & 
     220            IF( zcff  - nn_fsbc > 0.001_wp .OR. zcff  < 0._wp )   & 
    220221               &                   CALL ctl_stop( 'abl_init : wrong value for rn_ldyn_min' ) 
    221222            IF( zcff  > zcff1                    )   & 
    222                &                   CALL ctl_stop( 'abl_init : rn_ldyn_max should be smaller than rn_ldyn_min' ) 
     223               &                   CALL ctl_stop( 'abl_init : rn_ldyn_max must be smaller than rn_ldyn_min' ) 
    223224         END IF 
    224225      END IF 
     
    234235         ! Check that restoring coefficients are between 0 and 1 
    235236         !IF( zcff1 > 1._wp .OR. zcff1 < 0._wp )   & 
    236          IF( zcff1 > nn_fsbc .OR. zcff1 < 0._wp )   & 
     237         IF( zcff1 - nn_fsbc > 0.001_wp .OR. zcff1 < 0._wp )   & 
    237238            &                   CALL ctl_stop( 'abl_init : wrong value for rn_ltra_max' ) 
    238239         !IF( zcff  > 1._wp .OR. zcff  < 0._wp )   & 
    239          IF( zcff  > nn_fsbc .OR. zcff  < 0._wp )   & 
     240         IF( zcff  - nn_fsbc > 0.001_wp .OR. zcff  < 0._wp )   & 
    240241            &                   CALL ctl_stop( 'abl_init : wrong value for rn_ltra_min' ) 
    241242         IF( zcff  > zcff1                    )   & 
    242             &                   CALL ctl_stop( 'abl_init : rn_ltra_max should be smaller than rn_ltra_min' ) 
     243            &                   CALL ctl_stop( 'abl_init : rn_ltra_max must be smaller than rn_ltra_min' ) 
    243244      END IF 
    244245 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/ice.F90

    r11413 r11586  
    110110   !! bv_i        |      -      |    relative brine volume        | ???   |  
    111111   !! at_ip       |      -      |    Total ice pond concentration |       | 
     112   !! hm_ip       |      -      |    Mean ice pond depth          | m     | 
    112113   !! vt_ip       |      -      |    Total ice pond vol. per unit area| m | 
    113114   !!===================================================================== 
     
    188189 
    189190   !                                     !!** ice-ponds namelist (namthd_pnd) 
     191   LOGICAL , PUBLIC ::   ln_pnd           !: Melt ponds (T) or not (F) 
    190192   LOGICAL , PUBLIC ::   ln_pnd_H12       !: Melt ponds scheme from Holland et al 2012 
    191193   LOGICAL , PUBLIC ::   ln_pnd_CST       !: Melt ponds scheme with constant fraction and depth 
     
    196198   !                                     !!** ice-diagnostics namelist (namdia) ** 
    197199   LOGICAL , PUBLIC ::   ln_icediachk     !: flag for ice diag (T) or not (F) 
     200   REAL(wp), PUBLIC ::   rn_icechk_cel    !: rate of ice spuriously gained/lost (at any gridcell) 
     201   REAL(wp), PUBLIC ::   rn_icechk_glo    !: rate of ice spuriously gained/lost (globally) 
    198202   LOGICAL , PUBLIC ::   ln_icediahsb     !: flag for ice diag (T) or not (F) 
    199203   LOGICAL , PUBLIC ::   ln_icectl        !: flag for sea-ice points output (T) or not (F) 
     
    330334 
    331335   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   at_ip      !: total melt pond fraction 
     336   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hm_ip      !: mean melt pond depth                     [m] 
    332337   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   vt_ip      !: total melt pond volume per unit area     [m] 
    333338 
     
    351356   !! * Ice diagnostics 
    352357   !!---------------------------------------------------------------------- 
    353    ! thd refers to changes induced by thermodynamics 
    354    ! trp   ''         ''     ''       advection (transport of ice) 
    355    ! 
    356358   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_vi   !: transport of ice volume 
    357359   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_vs   !: transport of snw volume 
     
    365367   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_vsnw     !: snw volume variation   [m/s]  
    366368 
     369   !!---------------------------------------------------------------------- 
     370   !! * Ice conservation 
     371   !!---------------------------------------------------------------------- 
     372   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_v        !: conservation of ice volume 
     373   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_s        !: conservation of ice salt 
     374   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_t        !: conservation of ice heat 
     375   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_fv       !: conservation of ice volume 
     376   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_fs       !: conservation of ice salt 
     377   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_ft       !: conservation of ice heat 
    367378   ! 
    368379   !!---------------------------------------------------------------------- 
     
    389400      INTEGER :: ice_alloc 
    390401      ! 
    391       INTEGER :: ierr(15), ii 
     402      INTEGER :: ierr(16), ii 
    392403      !!----------------------------------------------------------------- 
    393404      ierr(:) = 0 
     
    440451 
    441452      ii = ii + 1 
    442       ALLOCATE( at_ip(jpi,jpj) , vt_ip(jpi,jpj) , STAT = ierr(ii) ) 
     453      ALLOCATE( at_ip(jpi,jpj) , hm_ip(jpi,jpj) , vt_ip(jpi,jpj) , STAT = ierr(ii) ) 
    443454 
    444455      ! * Old values of global variables 
     
    461472         &      diag_sice  (jpi,jpj) , diag_vice   (jpi,jpj) , diag_vsnw  (jpi,jpj), STAT=ierr(ii) ) 
    462473 
     474      ! * Ice conservation 
     475      ii = ii + 1 
     476      ALLOCATE( diag_v (jpi,jpj) , diag_s (jpi,jpj) , diag_t (jpi,jpj),   &  
     477         &      diag_fv(jpi,jpj) , diag_fs(jpi,jpj) , diag_ft(jpi,jpj), STAT=ierr(ii) ) 
     478       
    463479      ! * SIMIP diagnostics 
    464480      ii = ii + 1 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icecor.F90

    r11413 r11586  
    6060      IF( ln_timing    )   CALL timing_start('icecor')                                                             ! timing 
    6161      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icecor', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
     62      IF( ln_icediachk )   CALL ice_cons2D  (0, 'icecor',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation 
    6263      ! 
    6364      IF( kt == nit000 .AND. lwp .AND. kn == 2 ) THEN 
     
    164165      ! 
    165166      ! controls 
    166       IF( ln_icediachk   )   CALL ice_cons_hsm(1, 'icecor', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    167       IF( ln_ctl         )   CALL ice_prt3D   ('icecor')                                                             ! prints 
     167      IF( ln_ctl       )   CALL ice_prt3D   ('icecor')                                                             ! prints 
    168168      IF( ln_icectl .AND. kn == 2 ) & 
    169          &                   CALL ice_prt     ( kt, iiceprt, jiceprt, 2, ' - Final state - ' )                       ! prints 
    170       IF( ln_timing      )   CALL timing_stop ('icecor')                                                             ! timing 
     169         &                 CALL ice_prt     ( kt, iiceprt, jiceprt, 2, ' - Final state - ' )                       ! prints 
     170      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icecor', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
     171      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icecor',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation 
     172      IF( ln_timing    )   CALL timing_stop ('icecor')                                                             ! timing 
    171173      ! 
    172174   END SUBROUTINE ice_cor 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icectl.F90

    r11413 r11586  
    1212   !!   'key_si3'                                       SI3 sea-ice model 
    1313   !!---------------------------------------------------------------------- 
    14    !!    ice_cons_hsm     : conservation tests on heat, salt and mass 
    15    !!    ice_cons_final   : conservation tests on heat, salt and mass at end of time step 
     14   !!    ice_cons_hsm     : conservation tests on heat, salt and mass during a  time step (global)  
     15   !!    ice_cons_final   : conservation tests on heat, salt and mass at end of time step (global) 
     16   !!    ice_cons2D       : conservation tests on heat, salt and mass at each gridcell 
    1617   !!    ice_ctl          : control prints in case of crash 
    1718   !!    ice_prt          : control prints at a given grid point 
     
    2728   ! 
    2829   USE in_out_manager ! I/O manager 
     30   USE iom            ! I/O manager library 
    2931   USE lib_mpp        ! MPP library 
    3032   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero) 
     
    3739   PUBLIC   ice_cons_hsm 
    3840   PUBLIC   ice_cons_final 
     41   PUBLIC   ice_cons2D 
    3942   PUBLIC   ice_ctl 
    4043   PUBLIC   ice_prt 
    4144   PUBLIC   ice_prt3D 
    4245 
     46   ! thresold values for conservation 
     47   !    these values are changed by the namelist parameter rn_icechk, so that threshold = zchk * rn_icechk 
     48   REAL(wp), PARAMETER ::   zchk_m   = 1.e-5   ! kg/m2/s <=> 1mm of ice per year  spuriously gained/lost 
     49   REAL(wp), PARAMETER ::   zchk_s   = 1.e-4   ! g/m2/s  <=> 1mm of ice per year  spuriously gained/lost (considering s=10g/kg) 
     50   REAL(wp), PARAMETER ::   zchk_t   = 3.      ! W/m2    <=> 1mm of ice per year  spuriously gained/lost (considering Lf=3e5J/kg) 
     51    
    4352   !! * Substitutions 
    4453#  include "vectopt_loop_substitute.h90" 
     
    5968      !! ** Method  : This is an online diagnostics which can be activated with ln_icediachk=true 
    6069      !!              It prints in ocean.output if there is a violation of conservation at each time-step 
    61       !!              The thresholds (zv_sill, zs_sill, zt_sill) which determine violations are set to 
     70      !!              The thresholds (zchk_m, zchk_s, zchk_t) which determine violations are set to 
    6271      !!              a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years. 
    6372      !!              For salt and heat thresholds, ice is considered to have a salinity of 10  
     
    6877      REAL(wp)        , INTENT(inout) ::   pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft 
    6978      !! 
    70       REAL(wp) ::   zv, zs, zt, zfs, zfv, zft 
    71       REAL(wp) ::   zvmin, zamin, zamax, zeimin, zesmin, zsmin 
     79      REAL(wp) ::   zdiag_mass, zdiag_salt, zdiag_heat, & 
     80         &          zdiag_vmin, zdiag_amin, zdiag_amax, zdiag_eimin, zdiag_esmin, zdiag_smin 
    7281      REAL(wp) ::   zvtrp, zetrp 
    73       REAL(wp) ::   zarea, zv_sill, zs_sill, zt_sill 
    74       REAL(wp), PARAMETER ::   zconv = 1.e-9 ! convert W to GW and kg to Mt 
     82      REAL(wp) ::   zarea 
    7583      !!------------------------------------------------------------------- 
    7684      ! 
    7785      IF( icount == 0 ) THEN 
    78          !                          ! water flux 
    79          pdiag_fv = glob_sum( 'icectl',                                                                       & 
    80             &                 -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) +                  & 
    81             &                    wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_pnd(:,:)  +  & 
    82             &                    wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) +  & 
    83             &                    wfx_ice_sub(:,:) + wfx_spr(:,:)  & 
    84             &                  ) * e1e2t(:,:) ) * zconv 
     86 
     87         pdiag_v = glob_sum( 'icectl',   SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t ) 
     88         pdiag_s = glob_sum( 'icectl',   SUM( sv_i * rhoi            , dim=3 ) * e1e2t ) 
     89         pdiag_t = glob_sum( 'icectl', ( SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) ) * e1e2t ) 
     90 
     91         ! mass flux 
     92         pdiag_fv = glob_sum( 'icectl',  & 
     93            &                         ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd + & 
     94            &                           wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr ) * e1e2t ) 
     95         ! salt flux 
     96         pdiag_fs = glob_sum( 'icectl',  & 
     97            &                         ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + & 
     98            &                           sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) * e1e2t ) 
     99         ! heat flux 
     100         pdiag_ft = glob_sum( 'icectl',  & 
     101            &                         (   hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw  & 
     102            &                           - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) * e1e2t ) 
     103 
     104      ELSEIF( icount == 1 ) THEN 
     105 
     106         ! -- mass diag -- ! 
     107         zdiag_mass = ( glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t ) - pdiag_v ) * r1_rdtice       & 
     108            &         + glob_sum( 'icectl', ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn +       & 
     109            &                                 wfx_lam + wfx_pnd + wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + & 
     110            &                                 wfx_ice_sub + wfx_spr ) * e1e2t )                                           & 
     111            &         - pdiag_fv 
    85112         ! 
    86          !                          ! salt flux 
    87          pdiag_fs = glob_sum( 'icectl',                                                                     & 
    88             &                  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
    89             &                    sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)    & 
    90             &                  ) *  e1e2t(:,:) ) * zconv  
     113         ! -- salt diag -- ! 
     114         zdiag_salt = ( glob_sum( 'icectl', SUM( sv_i * rhoi , dim=3 ) * e1e2t ) - pdiag_s ) * r1_rdtice  & 
     115            &         + glob_sum( 'icectl', ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni +           & 
     116            &                                 sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) * e1e2t ) & 
     117            &         - pdiag_fs 
    91118         ! 
    92          !                          ! heat flux 
    93          pdiag_ft = glob_sum( 'icectl',                                                                    & 
    94             &                  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  &  
    95             &                  - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   & 
    96             &                  ) *  e1e2t(:,:) ) * zconv 
    97  
    98          pdiag_v = glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t * zconv ) 
    99  
    100          pdiag_s = glob_sum( 'icectl', SUM( sv_i * rhoi            , dim=3 ) * e1e2t * zconv ) 
    101  
    102          pdiag_t = glob_sum( 'icectl', (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )     & 
    103             &                 + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv 
    104  
    105       ELSEIF( icount == 1 ) THEN 
    106  
    107          ! water flux 
    108          zfv = glob_sum( 'icectl',                                                                        & 
    109             &             -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) +                  & 
    110             &                wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_pnd(:,:)  +  & 
    111             &                wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) +  & 
    112             &                wfx_ice_sub(:,:) + wfx_spr(:,:)  & 
    113             &              ) * e1e2t(:,:) ) * zconv - pdiag_fv 
    114  
    115          ! salt flux 
    116          zfs = glob_sum( 'icectl',                                                                       & 
    117             &              ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
    118             &                sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)    &  
    119             &              ) * e1e2t(:,:) ) * zconv - pdiag_fs 
    120  
    121          ! heat flux 
    122          zft = glob_sum( 'icectl',                                                                      & 
    123             &              ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  &  
    124             &              - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   & 
    125             &              ) * e1e2t(:,:) ) * zconv - pdiag_ft 
    126   
    127          ! outputs 
    128          zv = ( ( glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t ) * zconv  & 
    129             &     - pdiag_v ) * r1_rdtice - zfv ) * rday 
    130  
    131          zs = ( ( glob_sum( 'icectl', SUM( sv_i * rhoi            , dim=3 ) * e1e2t ) * zconv  & 
    132             &     - pdiag_s ) * r1_rdtice + zfs ) * rday 
    133  
    134          zt = ( glob_sum( 'icectl',                                                                & 
    135             &             (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )                       & 
    136             &              + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv   & 
    137             &   - pdiag_t ) * r1_rdtice + zft 
    138  
    139          ! zvtrp and zetrp must be close to 0 if the advection scheme is conservative 
    140          zvtrp = glob_sum( 'icectl', ( diag_trp_vi * rhoi + diag_trp_vs * rhos ) * e1e2t  ) * zconv * rday  
    141          zetrp = glob_sum( 'icectl', ( diag_trp_ei        + diag_trp_es        ) * e1e2t  ) * zconv 
    142  
    143          zamax  = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 
    144          zvmin  = glob_min( 'icectl', v_i ) 
    145          zamin  = glob_min( 'icectl', a_i ) 
    146          zsmin  = glob_min( 'icectl', sv_i ) 
    147          zeimin = glob_min( 'icectl', SUM( e_i, dim=3 ) ) 
    148          zesmin = glob_min( 'icectl', SUM( e_s, dim=3 ) ) 
    149  
    150          ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice)  
    151          zarea   = glob_sum( 'icectl', SUM( a_i + epsi10, dim=3 ) * e1e2t ) * zconv ! in 1.e9 m2 
    152          zv_sill = zarea * 2.5e-5 
    153          zs_sill = zarea * 25.e-5 
    154          zt_sill = zarea * 10.e-5 
    155  
    156          IF(lwp) THEN 
     119         ! -- heat diag -- ! 
     120         zdiag_heat = ( glob_sum( 'icectl', ( SUM(SUM(e_i, dim=4), dim=3) + SUM(SUM(e_s, dim=4), dim=3) ) * e1e2t ) - pdiag_t & 
     121            &         ) * r1_rdtice                                                                                           & 
     122            &         + glob_sum( 'icectl', (  hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw                      & 
     123            &                                - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) * e1e2t )                    & 
     124            &         - pdiag_ft 
     125 
     126         ! -- min/max diag -- ! 
     127         zdiag_amax  = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 
     128         zdiag_vmin  = glob_min( 'icectl', v_i ) 
     129         zdiag_amin  = glob_min( 'icectl', a_i ) 
     130         zdiag_smin  = glob_min( 'icectl', sv_i ) 
     131         zdiag_eimin = glob_min( 'icectl', SUM( e_i, dim=3 ) ) 
     132         zdiag_esmin = glob_min( 'icectl', SUM( e_s, dim=3 ) ) 
     133 
     134         ! -- advection scheme is conservative? -- ! 
     135         zvtrp = glob_sum( 'icectl', ( diag_trp_vi * rhoi + diag_trp_vs * rhos ) * e1e2t ) ! must be close to 0 
     136         zetrp = glob_sum( 'icectl', ( diag_trp_ei        + diag_trp_es        ) * e1e2t ) ! must be close to 0 
     137 
     138         ! ice area (+epsi10 to set a threshold > 0 when there is no ice)  
     139         zarea = glob_sum( 'icectl', SUM( a_i + epsi10, dim=3 ) * e1e2t ) 
     140 
     141         IF( lwp ) THEN 
    157142            ! check conservation issues 
    158             IF ( ABS( zv ) > zv_sill )   WRITE(numout,*) 'violation volume [Mt/day]     (',cd_routine,') = ',zv 
    159             IF ( ABS( zs ) > zs_sill )   WRITE(numout,*) 'violation saline [Mkg/day]    (',cd_routine,') = ',zs 
    160             IF ( ABS( zt ) > zt_sill )   WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',zt 
     143            IF( ABS(zdiag_mass) > zchk_m * rn_icechk_glo * zarea ) & 
     144               &                   WRITE(numout,*)   cd_routine,' : violation mass cons. [kg] = ',zdiag_mass * rdt_ice 
     145            IF( ABS(zdiag_salt) > zchk_s * rn_icechk_glo * zarea ) & 
     146               &                   WRITE(numout,*)   cd_routine,' : violation salt cons. [g]  = ',zdiag_salt * rdt_ice 
     147            IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zarea ) & 
     148               &                   WRITE(numout,*)   cd_routine,' : violation heat cons. [J]  = ',zdiag_heat * rdt_ice 
     149            ! check negative values 
     150            IF( zdiag_vmin  < 0. ) WRITE(numout,*)   cd_routine,' : violation v_i < 0         = ',zdiag_vmin 
     151            IF( zdiag_amin  < 0. ) WRITE(numout,*)   cd_routine,' : violation a_i < 0         = ',zdiag_amin 
     152            IF( zdiag_smin  < 0. ) WRITE(numout,*)   cd_routine,' : violation s_i < 0         = ',zdiag_smin 
     153            IF( zdiag_eimin < 0. ) WRITE(numout,*)   cd_routine,' : violation e_i < 0         = ',zdiag_eimin 
     154            IF( zdiag_esmin < 0. ) WRITE(numout,*)   cd_routine,' : violation e_s < 0         = ',zdiag_esmin 
    161155            ! check maximum ice concentration 
    162             IF ( zamax > MAX( rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' )  & 
    163                &                         WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax 
    164             ! check negative values 
    165             IF ( zvmin  < 0. )           WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',zvmin 
    166             IF ( zamin  < 0. )           WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin 
    167             IF ( zsmin  < 0. )           WRITE(numout,*) 'violation s_i<0               (',cd_routine,') = ',zsmin 
    168             IF ( zeimin < 0. )           WRITE(numout,*) 'violation e_i<0               (',cd_routine,') = ',zeimin 
    169             IF ( zesmin < 0. )           WRITE(numout,*) 'violation e_s<0               (',cd_routine,') = ',zesmin 
    170 !clem: the following check fails (I think...) 
    171 !            IF ( ABS(zvtrp ) > zv_sill .AND. cd_routine == 'icedyn_adv' ) THEN 
    172 !                                           WRITE(numout,*) 'violation vtrp [Mt/day]       (',cd_routine,') = ',zvtrp 
    173 !                                           WRITE(numout,*) 'violation etrp [GW]           (',cd_routine,') = ',zetrp 
    174 !            ENDIF 
     156            IF( zdiag_amax > MAX(rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' ) & 
     157               &                   WRITE(numout,*)   cd_routine,' : violation a_i > amax      = ',zdiag_amax 
     158            ! check if advection scheme is conservative 
     159            IF( ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 
     160               &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [kg] = ',zvtrp * rdt_ice 
    175161         ENDIF 
    176162         ! 
     
    179165   END SUBROUTINE ice_cons_hsm 
    180166 
    181  
    182167   SUBROUTINE ice_cons_final( cd_routine ) 
    183168      !!------------------------------------------------------------------- 
     
    188173      !! ** Method  : This is an online diagnostics which can be activated with ln_icediachk=true 
    189174      !!              It prints in ocean.output if there is a violation of conservation at each time-step 
    190       !!              The thresholds (zv_sill, zs_sill, zt_sill) which determine the violation are set to 
     175      !!              The thresholds (zchk_m, zchk_s, zchk_t) which determine the violation are set to 
    191176      !!              a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years. 
    192177      !!              For salt and heat thresholds, ice is considered to have a salinity of 10  
    193178      !!              and a heat content of 3e5 J/kg (=latent heat of fusion)  
    194179      !!------------------------------------------------------------------- 
    195       CHARACTER(len=*), INTENT(in)    :: cd_routine    ! name of the routine 
    196       REAL(wp)                        :: zqmass, zhfx, zsfx, zvfx 
    197       REAL(wp)                        :: zarea, zv_sill, zs_sill, zt_sill 
    198       REAL(wp), PARAMETER             :: zconv = 1.e-9 ! convert W to GW and kg to Mt 
     180      CHARACTER(len=*), INTENT(in) ::   cd_routine    ! name of the routine 
     181      REAL(wp) ::   zdiag_mass, zdiag_salt, zdiag_heat 
     182      REAL(wp) ::   zarea 
    199183      !!------------------------------------------------------------------- 
    200184 
    201185      ! water flux 
    202       zvfx  = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) * zconv * rday 
    203  
    204       ! salt flux 
    205       zsfx  = glob_sum( 'icectl', ( sfx + diag_sice ) * e1e2t ) * zconv * rday 
    206  
    207       ! heat flux 
     186      ! -- mass diag -- ! 
     187      zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) 
     188 
     189      ! -- salt diag -- ! 
     190      zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice ) * e1e2t ) 
     191 
     192      ! -- heat diag -- ! 
    208193      ! clem: not the good formulation 
    209       !!zhfx  = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat + hfx_thd + hfx_dyn + hfx_res + hfx_sub + hfx_spr  & 
    210       !!   &                        ) * e1e2t ) * zconv 
    211  
    212       ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice)  
    213       zarea   = glob_sum( 'icectl', SUM( a_i + epsi10, dim=3 ) * e1e2t ) * zconv ! in 1.e9 m2 
    214       zv_sill = zarea * 2.5e-5 
    215       zs_sill = zarea * 25.e-5 
    216       zt_sill = zarea * 10.e-5 
    217  
    218       IF(lwp) THEN 
    219          IF( ABS( zvfx ) > zv_sill )   WRITE(numout,*) 'violation vfx  [Mt/day]       (',cd_routine,') = ',zvfx 
    220          IF( ABS( zsfx ) > zs_sill )   WRITE(numout,*) 'violation sfx  [Mkg/day]      (',cd_routine,') = ',zsfx 
    221          !!IF( ABS( zhfx ) > zt_sill )   WRITE(numout,*) 'violation hfx  [GW]           (',cd_routine,') = ',zhfx 
     194      !!zdiag_heat  = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat + hfx_thd + hfx_dyn + hfx_res + hfx_sub + hfx_spr  & 
     195      !!   &                              ) * e1e2t ) 
     196 
     197      ! ice area (+epsi10 to set a threshold > 0 when there is no ice)  
     198      zarea = glob_sum( 'icectl', SUM( a_i + epsi10, dim=3 ) * e1e2t ) 
     199 
     200      IF( lwp ) THEN 
     201         IF( ABS(zdiag_mass) > zchk_m * rn_icechk_glo * zarea ) & 
     202            &                   WRITE(numout,*) cd_routine,' : violation mass cons. [kg] = ',zdiag_mass * rdt_ice 
     203         IF( ABS(zdiag_salt) > zchk_s * rn_icechk_glo * zarea ) & 
     204            &                   WRITE(numout,*) cd_routine,' : violation salt cons. [g]  = ',zdiag_salt * rdt_ice 
     205         !!IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zarea ) WRITE(numout,*) cd_routine,' : violation heat cons. [J]  = ',zdiag_heat * rdt_ice 
    222206      ENDIF 
    223207      ! 
    224208   END SUBROUTINE ice_cons_final 
    225209 
     210   SUBROUTINE ice_cons2D( icount, cd_routine, pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft ) 
     211      !!------------------------------------------------------------------- 
     212      !!                       ***  ROUTINE ice_cons2D *** 
     213      !! 
     214      !! ** Purpose : Test the conservation of heat, salt and mass for each ice routine 
     215      !!                     + test if ice concentration and volume are > 0 
     216      !! 
     217      !! ** Method  : This is an online diagnostics which can be activated with ln_icediachk=true 
     218      !!              It stops the code if there is a violation of conservation at any gridcell 
     219      !!------------------------------------------------------------------- 
     220      INTEGER         , INTENT(in) ::   icount        ! called at: =0 the begining of the routine, =1  the end 
     221      CHARACTER(len=*), INTENT(in) ::   cd_routine    ! name of the routine 
     222      REAL(wp)        , DIMENSION(jpi,jpj), INTENT(inout) ::   pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft 
     223      !! 
     224      REAL(wp), DIMENSION(jpi,jpj) ::   zdiag_mass, zdiag_salt, zdiag_heat, & 
     225         &                              zdiag_amin, zdiag_vmin, zdiag_smin, zdiag_emin !!, zdiag_amax   
     226      INTEGER ::   jl, jk 
     227      LOGICAL ::   ll_stop_m = .FALSE. 
     228      LOGICAL ::   ll_stop_s = .FALSE. 
     229      LOGICAL ::   ll_stop_t = .FALSE. 
     230      CHARACTER(len=120) ::   clnam   ! filename for the output 
     231      !!------------------------------------------------------------------- 
     232      ! 
     233      IF( icount == 0 ) THEN 
     234 
     235         pdiag_v = SUM( v_i  * rhoi + v_s * rhos, dim=3 ) 
     236         pdiag_s = SUM( sv_i * rhoi             , dim=3 ) 
     237         pdiag_t = SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) 
     238 
     239         ! mass flux 
     240         pdiag_fv = wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd  +  & 
     241            &       wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr 
     242         ! salt flux 
     243         pdiag_fs = sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam  
     244         ! heat flux 
     245         pdiag_ft =   hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw  &  
     246            &       - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr 
     247 
     248      ELSEIF( icount == 1 ) THEN 
     249 
     250         ! -- mass diag -- ! 
     251         zdiag_mass =   ( SUM( v_i * rhoi + v_s * rhos, dim=3 ) - pdiag_v ) * r1_rdtice                             & 
     252            &         + ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd + & 
     253            &             wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr )           & 
     254            &         - pdiag_fv 
     255         IF( MAXVAL( ABS(zdiag_mass) ) > zchk_m * rn_icechk_cel )   ll_stop_m = .TRUE. 
     256         ! 
     257         ! -- salt diag -- ! 
     258         zdiag_salt =   ( SUM( sv_i * rhoi , dim=3 ) - pdiag_s ) * r1_rdtice                                                  & 
     259            &         + ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) & 
     260            &         - pdiag_fs 
     261         IF( MAXVAL( ABS(zdiag_salt) ) > zchk_s * rn_icechk_cel )   ll_stop_s = .TRUE. 
     262         ! 
     263         ! -- heat diag -- ! 
     264         zdiag_heat =   ( SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) - pdiag_t ) * r1_rdtice & 
     265            &         + (  hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw                                &  
     266            &            - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr )                                        & 
     267            &         - pdiag_ft 
     268         IF( MAXVAL( ABS(zdiag_heat) ) > zchk_t * rn_icechk_cel )   ll_stop_t = .TRUE. 
     269         ! 
     270         ! -- other diags -- ! 
     271         ! a_i < 0 
     272         zdiag_amin(:,:) = 0._wp 
     273         DO jl = 1, jpl 
     274            WHERE( a_i(:,:,jl) < 0._wp )   zdiag_amin(:,:) = 1._wp 
     275         ENDDO 
     276         ! v_i < 0 
     277         zdiag_vmin(:,:) = 0._wp 
     278         DO jl = 1, jpl 
     279            WHERE( v_i(:,:,jl) < 0._wp )   zdiag_vmin(:,:) = 1._wp 
     280         ENDDO 
     281         ! s_i < 0 
     282         zdiag_smin(:,:) = 0._wp 
     283         DO jl = 1, jpl 
     284            WHERE( s_i(:,:,jl) < 0._wp )   zdiag_smin(:,:) = 1._wp 
     285         ENDDO 
     286         ! e_i < 0 
     287         zdiag_emin(:,:) = 0._wp 
     288         DO jl = 1, jpl 
     289            DO jk = 1, nlay_i 
     290               WHERE( e_i(:,:,jk,jl) < 0._wp )   zdiag_emin(:,:) = 1._wp 
     291            ENDDO 
     292         ENDDO 
     293         ! a_i > amax 
     294         !WHERE( SUM( a_i, dim=3 ) > ( MAX( rn_amax_n, rn_amax_s ) + epsi10 )   ;   zdiag_amax(:,:) = 1._wp 
     295         !ELSEWHERE                                                             ;   zdiag_amax(:,:) = 0._wp 
     296         !END WHERE 
     297 
     298         IF( ll_stop_m .OR. ll_stop_s .OR. ll_stop_t ) THEN 
     299            clnam = 'diag_ice_conservation_'//cd_routine 
     300            CALL ice_cons_wri( clnam, zdiag_mass, zdiag_salt, zdiag_heat, zdiag_amin, zdiag_vmin, zdiag_smin, zdiag_emin ) 
     301         ENDIF 
     302 
     303         IF( ll_stop_m )   CALL ctl_stop( 'STOP', cd_routine//': ice mass conservation issue' ) 
     304         IF( ll_stop_s )   CALL ctl_stop( 'STOP', cd_routine//': ice salt conservation issue' ) 
     305         IF( ll_stop_t )   CALL ctl_stop( 'STOP', cd_routine//': ice heat conservation issue' ) 
     306          
     307      ENDIF 
     308 
     309   END SUBROUTINE ice_cons2D 
     310 
     311   SUBROUTINE ice_cons_wri( cdfile_name, pdiag_mass, pdiag_salt, pdiag_heat, pdiag_amin, pdiag_vmin, pdiag_smin, pdiag_emin ) 
     312      !!--------------------------------------------------------------------- 
     313      !!                 ***  ROUTINE ice_cons_wri  *** 
     314      !!         
     315      !! ** Purpose :   create a NetCDF file named cdfile_name which contains  
     316      !!                the instantaneous fields when conservation issue occurs 
     317      !! 
     318      !! ** Method  :   NetCDF files using ioipsl 
     319      !!---------------------------------------------------------------------- 
     320      CHARACTER(len=*), INTENT( in ) ::   cdfile_name      ! name of the file created 
     321      REAL(wp), DIMENSION(:,:), INTENT( in ) ::   pdiag_mass, pdiag_salt, pdiag_heat, & 
     322         &                                        pdiag_amin, pdiag_vmin, pdiag_smin, pdiag_emin !!, pdiag_amax   
     323      !! 
     324      INTEGER ::   inum 
     325      !!---------------------------------------------------------------------- 
     326      !  
     327      IF(lwp) WRITE(numout,*) 
     328      IF(lwp) WRITE(numout,*) 'ice_cons_wri : single instantaneous ice state' 
     329      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~  named :', cdfile_name, '...nc' 
     330      IF(lwp) WRITE(numout,*)                 
     331 
     332      CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl ) 
     333       
     334      CALL iom_rstput( 0, 0, inum, 'cons_mass', pdiag_mass(:,:) , ktype = jp_r8 )    ! ice mass spurious lost/gain 
     335      CALL iom_rstput( 0, 0, inum, 'cons_salt', pdiag_salt(:,:) , ktype = jp_r8 )    ! ice salt spurious lost/gain 
     336      CALL iom_rstput( 0, 0, inum, 'cons_heat', pdiag_heat(:,:) , ktype = jp_r8 )    ! ice heat spurious lost/gain 
     337      ! other diags 
     338      CALL iom_rstput( 0, 0, inum, 'aneg_count', pdiag_amin(:,:) , ktype = jp_r8 )    !  
     339      CALL iom_rstput( 0, 0, inum, 'vneg_count', pdiag_vmin(:,:) , ktype = jp_r8 )    !  
     340      CALL iom_rstput( 0, 0, inum, 'sneg_count', pdiag_smin(:,:) , ktype = jp_r8 )    !  
     341      CALL iom_rstput( 0, 0, inum, 'eneg_count', pdiag_emin(:,:) , ktype = jp_r8 )    !  
     342       
     343      CALL iom_close( inum ) 
     344 
     345   END SUBROUTINE ice_cons_wri 
    226346    
    227347   SUBROUTINE ice_ctl( kt ) 
     
    246366      ialert_id = 2 ! reference number of this alert 
    247367      cl_alname(ialert_id) = ' Incompat vol and con         '    ! name of the alert 
    248  
    249368      DO jl = 1, jpl 
    250369         DO jj = 1, jpj 
    251370            DO ji = 1, jpi 
    252371               IF(  v_i(ji,jj,jl) /= 0._wp   .AND.   a_i(ji,jj,jl) == 0._wp   ) THEN 
    253                   !WRITE(numout,*) ' ALERTE 2 :   Incompatible volume and concentration ' 
    254                   !WRITE(numout,*) ' at_i     ', at_i(ji,jj) 
    255                   !WRITE(numout,*) ' Point - category', ji, jj, jl 
    256                   !WRITE(numout,*) ' a_i *** a_i_b   ', a_i      (ji,jj,jl), a_i_b  (ji,jj,jl) 
    257                   !WRITE(numout,*) ' v_i *** v_i_b   ', v_i      (ji,jj,jl), v_i_b  (ji,jj,jl) 
     372                  WRITE(numout,*) ' ALERTE 2 :   Incompatible volume and concentration ' 
    258373                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    259374               ENDIF 
     
    269384         DO ji = 1, jpi 
    270385            IF(   h_i(ji,jj,jl)  >  50._wp   ) THEN 
     386               WRITE(numout,*) ' ALERTE 3 :   Very thick ice' 
    271387               !CALL ice_prt( kt, ji, jj, 2, ' ALERTE 3 :   Very thick ice ' ) 
    272388               inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
     
    280396      DO jj = 1, jpj 
    281397         DO ji = 1, jpi 
    282             IF(   MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 1.5  .AND.  & 
     398            IF(   MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 2.  .AND.  & 
    283399               &  at_i(ji,jj) > 0._wp   ) THEN 
     400               WRITE(numout,*) ' ALERTE 4 :   Very fast ice' 
    284401               !CALL ice_prt( kt, ji, jj, 1, ' ALERTE 4 :   Very fast ice ' ) 
    285                !WRITE(numout,*) ' ice strength             : ', strength(ji,jj) 
    286                !WRITE(numout,*) ' oceanic stress utau      : ', utau(ji,jj)  
    287                !WRITE(numout,*) ' oceanic stress vtau      : ', vtau(ji,jj) 
    288                !WRITE(numout,*) ' sea-ice stress utau_ice  : ', utau_ice(ji,jj)  
    289                !WRITE(numout,*) ' sea-ice stress vtau_ice  : ', vtau_ice(ji,jj) 
    290                !WRITE(numout,*) ' sst                      : ', sst_m(ji,jj) 
    291                !WRITE(numout,*) ' sss                      : ', sss_m(ji,jj) 
    292                !WRITE(numout,*)  
     402               inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
     403            ENDIF 
     404         END DO 
     405      END DO 
     406 
     407      ! Alert on salt flux 
     408      ialert_id = 5 ! reference number of this alert 
     409      cl_alname(ialert_id) = ' High salt flux               ' ! name of the alert 
     410      DO jj = 1, jpj 
     411         DO ji = 1, jpi 
     412            IF( ABS( sfx (ji,jj) ) > 1.0e-2 ) THEN  ! = 1 psu/day for 1m ocean depth 
     413               WRITE(numout,*) ' ALERTE 5 :   High salt flux' 
     414               !CALL ice_prt( kt, ji, jj, 3, ' ALERTE 5 :   High salt flux ' ) 
    293415               inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    294416            ENDIF 
     
    302424         DO ji = 1, jpi 
    303425            IF(   tmask(ji,jj,1) <= 0._wp   .AND.   at_i(ji,jj) > 0._wp   ) THEN  
     426               WRITE(numout,*) ' ALERTE 6 :   Ice on continents' 
    304427               !CALL ice_prt( kt, ji, jj, 1, ' ALERTE 6 :   Ice on continents ' ) 
    305                !WRITE(numout,*) ' masks s, u, v        : ', tmask(ji,jj,1), umask(ji,jj,1), vmask(ji,jj,1)  
    306                !WRITE(numout,*) ' sst                  : ', sst_m(ji,jj) 
    307                !WRITE(numout,*) ' sss                  : ', sss_m(ji,jj) 
    308                !WRITE(numout,*) ' at_i(ji,jj)          : ', at_i(ji,jj) 
    309                !WRITE(numout,*) ' v_ice(ji,jj)         : ', v_ice(ji,jj) 
    310                !WRITE(numout,*) ' v_ice(ji,jj-1)       : ', v_ice(ji,jj-1) 
    311                !WRITE(numout,*) ' u_ice(ji-1,jj)       : ', u_ice(ji-1,jj) 
    312                !WRITE(numout,*) ' u_ice(ji,jj)         : ', v_ice(ji,jj) 
    313                ! 
    314428               inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    315429            ENDIF 
     
    325439            DO ji = 1, jpi 
    326440               IF( s_i(ji,jj,jl) < 0.1 .AND. a_i(ji,jj,jl) > 0._wp ) THEN 
     441                  WRITE(numout,*) ' ALERTE 7 :   Very fresh ice' 
    327442!                 CALL ice_prt(kt,ji,jj,1, ' ALERTE 7 :   Very fresh ice ' ) 
    328 !                 WRITE(numout,*) ' sst                  : ', sst_m(ji,jj) 
    329 !                 WRITE(numout,*) ' sss                  : ', sss_m(ji,jj) 
    330 !                 WRITE(numout,*)  
    331443                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    332444               ENDIF 
     
    335447      END DO 
    336448! 
     449      ! Alert if qns very big 
     450      ialert_id = 8 ! reference number of this alert 
     451      cl_alname(ialert_id) = ' fnsolar very big             ' ! name of the alert 
     452      DO jj = 1, jpj 
     453         DO ji = 1, jpi 
     454            IF( ABS( qns(ji,jj) ) > 1500._wp .AND. at_i(ji,jj) > 0._wp ) THEN 
     455               ! 
     456               WRITE(numout,*) ' ALERTE 8 :   Very high non-solar heat flux' 
     457               !CALL ice_prt( kt, ji, jj, 2, '   ') 
     458               inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
     459               ! 
     460            ENDIF 
     461         END DO 
     462      END DO 
     463      !+++++ 
    337464 
    338465!     ! Alert if too old ice 
     
    345472                      ( ABS( o_i(ji,jj,jl) ) < 0._wp) ) .AND. & 
    346473                             ( a_i(ji,jj,jl) > 0._wp ) ) THEN 
     474                  WRITE(numout,*) ' ALERTE 9 :   Wrong ice age' 
    347475                  !CALL ice_prt( kt, ji, jj, 1, ' ALERTE 9 :   Wrong ice age ') 
    348476                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
     
    351479         END DO 
    352480      END DO 
    353   
    354       ! Alert on salt flux 
    355       ialert_id = 5 ! reference number of this alert 
    356       cl_alname(ialert_id) = ' High salt flux               ' ! name of the alert 
    357       DO jj = 1, jpj 
    358          DO ji = 1, jpi 
    359             IF( ABS( sfx (ji,jj) ) > 1.0e-2 ) THEN  ! = 1 psu/day for 1m ocean depth 
    360                !CALL ice_prt( kt, ji, jj, 3, ' ALERTE 5 :   High salt flux ' ) 
    361                !DO jl = 1, jpl 
    362                   !WRITE(numout,*) ' Category no: ', jl 
    363                   !WRITE(numout,*) ' a_i        : ', a_i      (ji,jj,jl) , ' a_i_b      : ', a_i_b  (ji,jj,jl)    
    364                   !WRITE(numout,*) ' v_i        : ', v_i      (ji,jj,jl) , ' v_i_b      : ', v_i_b  (ji,jj,jl)    
    365                   !WRITE(numout,*) ' ' 
    366                !END DO 
    367                inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    368             ENDIF 
    369          END DO 
    370       END DO 
    371  
    372       ! Alert if qns very big 
    373       ialert_id = 8 ! reference number of this alert 
    374       cl_alname(ialert_id) = ' fnsolar very big             ' ! name of the alert 
    375       DO jj = 1, jpj 
    376          DO ji = 1, jpi 
    377             IF( ABS( qns(ji,jj) ) > 1500._wp .AND. at_i(ji,jj) > 0._wp ) THEN 
    378                ! 
    379                !WRITE(numout,*) ' ALERTE 8 :   Very high non-solar heat flux' 
    380                !WRITE(numout,*) ' ji, jj    : ', ji, jj 
    381                !WRITE(numout,*) ' qns       : ', qns(ji,jj) 
    382                !WRITE(numout,*) ' sst       : ', sst_m(ji,jj) 
    383                !WRITE(numout,*) ' sss       : ', sss_m(ji,jj) 
    384                ! 
    385                !CALL ice_prt( kt, ji, jj, 2, '   ') 
    386                inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    387                ! 
    388             ENDIF 
    389          END DO 
    390       END DO 
    391       !+++++ 
    392   
     481   
    393482      ! Alert if very warm ice 
    394483      ialert_id = 10 ! reference number of this alert 
     
    400489               DO ji = 1, jpi 
    401490                  ztmelts    =  -rTmlt * sz_i(ji,jj,jk,jl) + rt0 
    402                   IF( t_i(ji,jj,jk,jl) >= ztmelts  .AND.  v_i(ji,jj,jl) > 1.e-10   & 
    403                      &                             .AND.  a_i(ji,jj,jl) > 0._wp   ) THEN 
    404                      !WRITE(numout,*) ' ALERTE 10 :   Very warm ice' 
    405                      !WRITE(numout,*) ' ji, jj, jk, jl : ', ji, jj, jk, jl 
    406                      !WRITE(numout,*) ' t_i : ', t_i(ji,jj,jk,jl) 
    407                      !WRITE(numout,*) ' e_i : ', e_i(ji,jj,jk,jl) 
    408                      !WRITE(numout,*) ' sz_i: ', sz_i(ji,jj,jk,jl) 
    409                      !WRITE(numout,*) ' ztmelts : ', ztmelts 
    410                      inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
     491                  IF( t_i(ji,jj,jk,jl) > ztmelts  .AND.  v_i(ji,jj,jl) > 1.e-10   & 
     492                     &                            .AND.  a_i(ji,jj,jl) > 0._wp   ) THEN 
     493                     WRITE(numout,*) ' ALERTE 10 :   Very warm ice' 
     494                    inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    411495                  ENDIF 
    412496               END DO 
     
    435519   END SUBROUTINE ice_ctl 
    436520  
    437     
    438521   SUBROUTINE ice_prt( kt, ki, kj, kn, cd1 ) 
    439522      !!------------------------------------------------------------------- 
     
    443526      !!                in ocean.ouput  
    444527      !!                3 possibilities exist  
    445       !!                n = 1/-1 -> simple ice state (plus Mechanical Check if -1) 
     528      !!                n = 1/-1 -> simple ice state 
    446529      !!                n = 2    -> exhaustive state 
    447530      !!                n = 3    -> ice/ocean salt fluxes 
     
    482565               WRITE(numout,*) ' - Cell values ' 
    483566               WRITE(numout,*) '   ~~~~~~~~~~~ ' 
    484                WRITE(numout,*) ' cell area     : ', e1e2t(ji,jj) 
    485567               WRITE(numout,*) ' at_i          : ', at_i(ji,jj)        
     568               WRITE(numout,*) ' ato_i         : ', ato_i(ji,jj)        
    486569               WRITE(numout,*) ' vt_i          : ', vt_i(ji,jj)        
    487570               WRITE(numout,*) ' vt_s          : ', vt_s(ji,jj)        
     
    503586               END DO 
    504587            ENDIF 
    505             IF( kn == -1 ) THEN 
    506                WRITE(numout,*) ' Mechanical Check ************** ' 
    507                WRITE(numout,*) ' Check what means ice divergence ' 
    508                WRITE(numout,*) ' Total ice concentration ', at_i (ji,jj) 
    509                WRITE(numout,*) ' Total lead fraction     ', ato_i(ji,jj) 
    510                WRITE(numout,*) ' Sum of both             ', ato_i(ji,jj) + at_i(ji,jj) 
    511                WRITE(numout,*) ' Sum of both minus 1     ', ato_i(ji,jj) + at_i(ji,jj) - 1.00 
    512             ENDIF 
    513              
    514588 
    515589            !-------------------- 
     
    525599               WRITE(numout,*) ' - Cell values ' 
    526600               WRITE(numout,*) '   ~~~~~~~~~~~ ' 
    527                WRITE(numout,*) ' cell area     : ', e1e2t(ji,jj) 
    528601               WRITE(numout,*) ' at_i          : ', at_i(ji,jj)        
    529602               WRITE(numout,*) ' vt_i          : ', vt_i(ji,jj)        
     
    624697      !! 
    625698      !!------------------------------------------------------------------- 
    626       CHARACTER(len=*), INTENT(in)  :: cd_routine    ! name of the routine 
    627       INTEGER                       :: jk, jl        ! dummy loop indices 
     699      CHARACTER(len=*), INTENT(in) ::  cd_routine    ! name of the routine 
     700      INTEGER                      ::  jk, jl        ! dummy loop indices 
    628701       
    629702      CALL prt_ctl_info(' ========== ') 
     
    684757       
    685758   END SUBROUTINE ice_prt3D 
    686  
     759       
    687760#else 
    688761   !!---------------------------------------------------------------------- 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icedia.F90

    r11413 r11586  
    175175      INTEGER            ::   ios, ierror   ! local integer 
    176176      !! 
    177       NAMELIST/namdia/ ln_icediachk, ln_icediahsb, ln_icectl, iiceprt, jiceprt   
     177      NAMELIST/namdia/ ln_icediachk, rn_icechk_cel, rn_icechk_glo, ln_icediahsb, ln_icectl, iiceprt, jiceprt   
    178178      !!---------------------------------------------------------------------- 
    179179      ! 
     
    191191         WRITE(numout,*) ' ~~~~~~~~~~~' 
    192192         WRITE(numout,*) '   Namelist namdia:' 
    193          WRITE(numout,*) '      Diagnose online heat/mass/salt budget      ln_icediachk = ', ln_icediachk 
    194          WRITE(numout,*) '      Output          heat/mass/salt budget      ln_icediahsb = ', ln_icediahsb 
    195          WRITE(numout,*) '      control prints for a given grid point      ln_icectl    = ', ln_icectl 
    196          WRITE(numout,*) '         chosen grid point position         (iiceprt,jiceprt) = (', iiceprt,',', jiceprt,')' 
     193         WRITE(numout,*) '      Diagnose online heat/mass/salt conservation ln_icediachk  = ', ln_icediachk 
     194         WRITE(numout,*) '         threshold for conservation (gridcell)    rn_icechk_cel = ', rn_icechk_cel 
     195         WRITE(numout,*) '         threshold for conservation (global)      rn_icechk_glo = ', rn_icechk_glo 
     196         WRITE(numout,*) '      Output          heat/mass/salt budget       ln_icediahsb  = ', ln_icediahsb 
     197         WRITE(numout,*) '      control prints for a given grid point       ln_icectl     = ', ln_icectl 
     198         WRITE(numout,*) '         chosen grid point position          (iiceprt,jiceprt)  = (', iiceprt,',', jiceprt,')' 
    197199      ENDIF 
    198200      !       
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icedyn_rdgrft.F90

    r11348 r11586  
    145145      IF( ln_timing    )   CALL timing_start('icedyn_rdgrft')                                                             ! timing 
    146146      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
     147      IF( ln_icediachk )   CALL ice_cons2D  (0, 'icedyn_rdgrft',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation 
    147148 
    148149      IF( kt == nit000 ) THEN 
     
    276277 
    277278      ! controls 
     279      IF( ln_ctl       )   CALL ice_prt3D   ('icedyn_rdgrft')                                                             ! prints 
    278280      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    279       IF( ln_ctl       )   CALL ice_prt3D   ('icedyn_rdgrft')                                                             ! prints 
     281      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icedyn_rdgrft',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation 
    280282      IF( ln_timing    )   CALL timing_stop ('icedyn_rdgrft')                                                             ! timing 
    281283      ! 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icedyn_rhg.F90

    r11413 r11586  
    6464      IF( ln_timing    )   CALL timing_start('icedyn_rhg')                                                             ! timing 
    6565      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn_rhg', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
     66      IF( ln_icediachk )   CALL ice_cons2D  (0, 'icedyn_rhg',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation 
    6667      ! 
    6768      IF( kt == nit000 .AND. lwp ) THEN 
     
    8788      ! 
    8889      ! controls 
     90      IF( ln_ctl       )   CALL ice_prt3D   ('icedyn_rhg')                                                             ! prints 
    8991      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icedyn_rhg', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    90       IF( ln_ctl       )   CALL ice_prt3D   ('icedyn_rhg')                                                             ! prints 
     92      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icedyn_rhg',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation 
    9193      IF( ln_timing    )   CALL timing_stop ('icedyn_rhg')                                                             ! timing 
    9294      ! 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icedyn_rhg_evp.F90

    r11413 r11586  
    121121      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV, zvU, zvV             ! ice/snow mass and volume 
    122122      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2       ! temporary scalars 
    123       REAL(wp) ::   zTauO, zTauB, zTauE, zvel                           ! temporary scalars 
     123      REAL(wp) ::   zTauO, zTauB, zRHS, zvel                            ! temporary scalars 
    124124      REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice 
    125125      REAL(wp) ::   zvCr                                                ! critical ice volume above which ice is landfast 
     
    130130      REAL(wp) ::   zshear, zdum1, zdum2 
    131131      ! 
    132       REAL(wp), DIMENSION(jpi,jpj) ::   z1_e1t0, z1_e2t0                ! scale factors 
    133132      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points 
    134133      REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017 
     
    244243      CALL ice_strength 
    245244 
    246       ! scale factors 
    247       DO jj = 2, jpjm1 
    248          DO ji = fs_2, fs_jpim1 
    249             z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj  ) + e1t(ji,jj  ) ) 
    250             z1_e2t0(ji,jj) = 1._wp / ( e2t(ji  ,jj+1) + e2t(ji,jj  ) ) 
    251          END DO 
    252       END DO 
    253  
    254245      ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 
    255246      IF( ln_landfast_L16 ) THEN   ;   zkt = rn_tensile 
     
    280271 
    281272            ! Ocean currents at U-V points 
    282             v_oceU(ji,jj)   = 0.5_wp * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji+1,jj)    & 
    283                &                       + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
    284              
    285             u_oceV(ji,jj)   = 0.5_wp * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj+1)    & 
    286                &                       + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
     273            v_oceU(ji,jj)   = 0.25_wp * ( v_oce(ji,jj) + v_oce(ji,jj-1) + v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * umask(ji,jj,1) 
     274            u_oceV(ji,jj)   = 0.25_wp * ( u_oce(ji,jj) + u_oce(ji-1,jj) + u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * vmask(ji,jj,1) 
    287275 
    288276            ! Coriolis at T points (m*f) 
     
    459447                  &                  ) * r1_e1e2v(ji,jj) 
    460448               ! 
    461                !                !--- u_ice at V point 
    462                u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
    463                   &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
     449               !                !--- ice currents at U-V point 
     450               v_iceU(ji,jj) = 0.25_wp * ( v_ice(ji,jj) + v_ice(ji,jj-1) + v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * umask(ji,jj,1) 
     451               u_iceV(ji,jj) = 0.25_wp * ( u_ice(ji,jj) + u_ice(ji-1,jj) + u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * vmask(ji,jj,1) 
    464452               ! 
    465                !                !--- v_ice at U point 
    466                v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     & 
    467                   &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
    468453            END DO 
    469454         END DO 
     
    494479                  ! 
    495480                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    496                   zTauE = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
    497                   ! 
    498                   !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    499                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztauy_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     481                  zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
     482                  ! 
     483                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS) 
     484                  !                                         1 = sliding friction : TauB < RHS 
     485                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 
    500486                  ! 
    501487                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    502488                     v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity 
    503                         &                                  + zTauE + zTauO * v_ice(ji,jj) )                                        & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     489                        &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    504490                        &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    505491                        &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
     
    508494                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    509495                     v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity 
    510                         &                                     + zTauE + zTauO * v_ice(ji,jj) )                                     & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     496                        &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    511497                        &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
    512498                        &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
     
    544530                  ! 
    545531                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    546                   zTauE = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
    547                   ! 
    548                   !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    549                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztaux_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     532                  zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
     533                  ! 
     534                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS) 
     535                  !                                         1 = sliding friction : TauB < RHS 
     536                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 
    550537                  ! 
    551538                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    552539                     u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity 
    553                         &                                  + zTauE + zTauO * u_ice(ji,jj) )                                        & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     540                        &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    554541                        &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    555                         &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
    556                         &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    557                         &            )   * zmsk00x(ji,jj) 
     542                        &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
     543                        &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
     544                        &           )   * zmsk00x(ji,jj) 
    558545                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    559546                     u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity 
    560                         &                                     + zTauE + zTauO * u_ice(ji,jj) )                                     & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     547                        &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    561548                        &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
    562549                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
     
    596583                  ! 
    597584                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    598                   zTauE = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
    599                   ! 
    600                   !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    601                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztaux_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     585                  zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
     586                  ! 
     587                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS) 
     588                  !                                         1 = sliding friction : TauB < RHS 
     589                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 
    602590                  ! 
    603591                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    604592                     u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity 
    605                         &                                  + zTauE + zTauO * u_ice(ji,jj) )                                        & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     593                        &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    606594                        &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    607                         &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
    608                         &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    609                         &            )   * zmsk00x(ji,jj) 
     595                        &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
     596                        &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
     597                        &           )   * zmsk00x(ji,jj) 
    610598                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    611599                     u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity 
    612                         &                                     + zTauE + zTauO * u_ice(ji,jj) )                                     & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     600                        &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    613601                        &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
    614602                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
     
    646634                  ! 
    647635                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    648                   zTauE = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
    649                   ! 
    650                   !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    651                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE + ztauy_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     636                  zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
     637                  ! 
     638                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS) 
     639                  !                                         1 = sliding friction : TauB < RHS 
     640                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 
    652641                  ! 
    653642                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    654643                     v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity 
    655                         &                                  + zTauE + zTauO * v_ice(ji,jj) )                                        & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     644                        &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    656645                        &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    657646                        &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
     
    660649                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    661650                     v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity 
    662                         &                                     + zTauE + zTauO * v_ice(ji,jj) )                                     & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     651                        &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    663652                        &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
    664653                        &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/iceistate.F90

    r11413 r11586  
    101101      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zti_3d , zts_3d                            !temporary arrays 
    102102      !! 
    103       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d 
     103      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d 
    104104      !-------------------------------------------------------------------- 
    105105 
     
    209209            ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1) 
    210210            ! 
    211             ! ponds 
     211            ! pond concentration 
    212212            IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) & 
    213                &     si(jp_apd)%fnow(:,:,1) = ( rn_apd_ini_n * zswitch + rn_apd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 
     213               &     si(jp_apd)%fnow(:,:,1) = ( rn_apd_ini_n * zswitch + rn_apd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) & ! rn_apd = pond fraction => rn_apnd * a_i = pond conc. 
     214               &                              * si(jp_ati)%fnow(:,:,1)  
    214215            zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1) 
     216            ! 
     217            ! pond depth 
    215218            IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) & 
    216219               &     si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 
     
    238241               zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:) 
    239242               ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:) 
    240                zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) 
     243               zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.  
    241244               zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:) 
    242245            ELSEWHERE 
     
    248251               zt_su_ini(:,:) = rn_tsu_ini_s * zswitch(:,:) 
    249252               ztm_s_ini(:,:) = rn_tms_ini_s * zswitch(:,:) 
    250                zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) 
     253               zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 
    251254               zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:) 
    252255            END WHERE 
    253256            ! 
    254257         ENDIF 
     258 
     259         ! make sure ponds = 0 if no ponds scheme 
     260         IF ( .NOT.ln_pnd ) THEN 
     261            zapnd_ini(:,:) = 0._wp 
     262            zhpnd_ini(:,:) = 0._wp 
     263         ENDIF 
     264          
    255265         !-------------! 
    256266         ! fill fields ! 
     
    275285         CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti)  , zt_su_ini ) 
    276286         CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti)  , zsm_i_ini ) 
     287         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti)  , zapnd_ini ) 
     288         CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d(1:npti)  , zhpnd_ini ) 
    277289 
    278290         ! allocate temporary arrays 
    279291         ALLOCATE( zhi_2d(npti,jpl), zhs_2d(npti,jpl), zai_2d (npti,jpl), & 
    280             &      zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl) ) 
     292            &      zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), zaip_2d(npti,jpl), zhip_2d(npti,jpl) ) 
    281293          
    282294         ! distribute 1-cat into jpl-cat: (jpi*jpj) -> (jpi*jpj,jpl) 
    283          CALL ice_var_itd( h_i_1d(1:npti)  , h_s_1d(1:npti)  , at_i_1d(1:npti),                   zhi_2d, zhs_2d, zai_2d , & 
    284             &              t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti), s_i_1d(1:npti),   zti_2d, zts_2d, ztsu_2d, zsi_2d ) 
     295         CALL ice_var_itd( h_i_1d(1:npti)  , h_s_1d(1:npti)  , at_i_1d(1:npti),                                                   & 
     296            &              zhi_2d          , zhs_2d          , zai_2d         ,                                                   & 
     297            &              t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti), s_i_1d(1:npti), a_ip_1d(1:npti), h_ip_1d(1:npti), & 
     298            &              zti_2d          , zts_2d          , ztsu_2d        , zsi_2d        , zaip_2d        , zhip_2d ) 
    285299 
    286300         ! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl) 
     
    289303            zts_3d(:,:,jl) = rt0 * tmask(:,:,1) 
    290304         END DO 
    291          CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d  , h_i  ) 
    292          CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d  , h_s  ) 
    293          CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d  , a_i  ) 
    294          CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d  , zti_3d ) 
    295          CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d  , zts_3d ) 
    296          CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d , t_su ) 
    297          CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d  , s_i  ) 
     305         CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d   , h_i    ) 
     306         CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d   , h_s    ) 
     307         CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d   , a_i    ) 
     308         CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d   , zti_3d ) 
     309         CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d   , zts_3d ) 
     310         CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d  , t_su   ) 
     311         CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d   , s_i    ) 
     312         CALL tab_2d_3d( npti, nptidx(1:npti), zaip_2d  , a_ip   ) 
     313         CALL tab_2d_3d( npti, nptidx(1:npti), zhip_2d  , h_ip   ) 
    298314 
    299315         ! deallocate temporary arrays 
    300316         DEALLOCATE( zhi_2d, zhs_2d, zai_2d , & 
    301             &        zti_2d, zts_2d, ztsu_2d, zsi_2d ) 
    302  
    303          ! Melt ponds: distribute uniformely over the categories 
    304          IF ( ln_pnd_CST .OR. ln_pnd_H12 ) THEN 
    305             DO jl = 1, jpl 
    306                a_ip_frac(:,:,jl) = zapnd_ini(:,:) 
    307                h_ip     (:,:,jl) = zhpnd_ini(:,:) 
    308                a_ip     (:,:,jl) = a_ip_frac(:,:,jl) * a_i (:,:,jl)  
    309                v_ip     (:,:,jl) = h_ip     (:,:,jl) * a_ip(:,:,jl) 
    310             END DO 
    311          ENDIF 
    312            
     317            &        zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d ) 
     318 
    313319         ! calculate extensive and intensive variables 
    314320         CALL ice_var_salprof ! for sz_i 
     
    350356         END DO 
    351357 
     358         ! Melt ponds 
     359         WHERE( a_i > epsi10 ) 
     360            a_ip_frac(:,:,:) = a_ip(:,:,:) / a_i(:,:,:) 
     361         ELSEWHERE 
     362            a_ip_frac(:,:,:) = 0._wp 
     363         END WHERE 
     364         v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:) 
     365           
    352366         ! specific temperatures for coupled runs 
    353367         tn_ice(:,:,:) = t_su(:,:,:) 
     
    512526      ENDIF 
    513527      ! 
     528      IF( .NOT.ln_pnd ) THEN 
     529         rn_apd_ini_n = 0. ; rn_apd_ini_s = 0. 
     530         rn_hpd_ini_n = 0. ; rn_hpd_ini_s = 0. 
     531         CALL ctl_warn( 'rn_apd_ini & rn_hpd_ini = 0 when no ponds' ) 
     532      ENDIF 
     533      ! 
    514534   END SUBROUTINE ice_istate_init 
    515535 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/iceitd.F90

    r11348 r11586  
    8888 
    8989      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
     90      IF( ln_icediachk )   CALL ice_cons2D  (0, 'iceitd_rem',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) 
    9091 
    9192      !----------------------------------------------------------------------------------------------- 
     
    316317      ! 
    317318      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
     319      IF( ln_icediachk )   CALL ice_cons2D  (1, 'iceitd_rem',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) 
    318320      ! 
    319321   END SUBROUTINE ice_itd_rem 
     
    586588      ! 
    587589      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
     590      IF( ln_icediachk )   CALL ice_cons2D  (0, 'iceitd_reb',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) 
    588591      ! 
    589592      jdonor(:,:) = 0 
     
    664667      ! 
    665668      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
     669      IF( ln_icediachk )   CALL ice_cons2D  (1, 'iceitd_reb',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) 
    666670      ! 
    667671   END SUBROUTINE ice_itd_reb 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icethd.F90

    r11348 r11586  
    9595      IF( ln_timing    )   CALL timing_start('icethd')                                                             ! timing 
    9696      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
     97      IF( ln_icediachk )   CALL ice_cons2D  (0, 'icethd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation 
    9798 
    9899      IF( kt == nit000 .AND. lwp ) THEN 
     
    243244      ! 
    244245      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
     246      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icethd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) 
    245247      !                    
    246248      IF( jpl > 1  )          CALL ice_itd_rem( kt )                ! --- Transport ice between thickness categories --- ! 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icethd_do.F90

    r11348 r11586  
    113113 
    114114      IF( ln_icediachk )   CALL ice_cons_hsm( 0, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft ) 
     115      IF( ln_icediachk )   CALL ice_cons2D  ( 0, 'icethd_do',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft ) 
    115116 
    116117      at_i(:,:) = SUM( a_i, dim=3 ) 
     
    420421      ! 
    421422      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
     423      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icethd_do',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) 
    422424      ! 
    423425   END SUBROUTINE ice_thd_do 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icethd_pnd.F90

    r11348 r11586  
    205205      INTEGER  ::   ios, ioptio   ! Local integer 
    206206      !! 
    207       NAMELIST/namthd_pnd/  ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb 
     207      NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb 
    208208      !!------------------------------------------------------------------- 
    209209      ! 
     
    221221         WRITE(numout,*) '~~~~~~~~~~~~~~~~' 
    222222         WRITE(numout,*) '   Namelist namicethd_pnd:' 
    223          WRITE(numout,*) '      Evolutive  melt pond fraction and depth (Holland et al 2012) ln_pnd_H12 = ', ln_pnd_H12 
    224          WRITE(numout,*) '      Prescribed melt pond fraction and depth                      ln_pnd_CST = ', ln_pnd_CST 
    225          WRITE(numout,*) '         Prescribed pond fraction                                  rn_apnd    = ', rn_apnd 
    226          WRITE(numout,*) '         Prescribed pond depth                                     rn_hpnd    = ', rn_hpnd 
    227          WRITE(numout,*) '      Melt ponds affect albedo or not                              ln_pnd_alb = ', ln_pnd_alb 
     223         WRITE(numout,*) '      Melt ponds activated or not                                     ln_pnd     = ', ln_pnd 
     224         WRITE(numout,*) '         Evolutive  melt pond fraction and depth (Holland et al 2012) ln_pnd_H12 = ', ln_pnd_H12 
     225         WRITE(numout,*) '         Prescribed melt pond fraction and depth                      ln_pnd_CST = ', ln_pnd_CST 
     226         WRITE(numout,*) '            Prescribed pond fraction                                  rn_apnd    = ', rn_apnd 
     227         WRITE(numout,*) '            Prescribed pond depth                                     rn_hpnd    = ', rn_hpnd 
     228         WRITE(numout,*) '         Melt ponds affect albedo or not                              ln_pnd_alb = ', ln_pnd_alb 
    228229      ENDIF 
    229230      ! 
    230231      !                             !== set the choice of ice pond scheme ==! 
    231232      ioptio = 0 
    232                                                             nice_pnd = np_pndNO 
    233       IF( ln_pnd_CST ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF 
    234       IF( ln_pnd_H12 ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndH12    ;   ENDIF 
    235       IF( ioptio > 1 )   CALL ctl_stop( 'ice_thd_pnd_init: choose one and only one pond scheme (ln_pnd_H12 or ln_pnd_CST)' ) 
     233      IF( .NOT.ln_pnd ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndNO     ;   ENDIF 
     234      IF( ln_pnd_CST  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF 
     235      IF( ln_pnd_H12  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndH12    ;   ENDIF 
     236      IF( ioptio /= 1 )   & 
     237         & CALL ctl_stop( 'ice_thd_pnd_init: choose either none (ln_pnd=F) or only one pond scheme (ln_pnd_H12 or ln_pnd_CST)' ) 
    236238      ! 
    237239      SELECT CASE( nice_pnd ) 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icevar.F90

    r11413 r11586  
    4747   !!   ice_var_zapneg    : remove negative ice fields 
    4848   !!   ice_var_roundoff  : remove negative values arising from roundoff erros 
    49    !!   ice_var_itd       : convert N-cat to M-cat 
    5049   !!   ice_var_bv        : brine volume 
    5150   !!   ice_var_enthalpy  : compute ice and snow enthalpies from temperature 
    5251   !!   ice_var_sshdyn    : compute equivalent ssh in lead 
     52   !!   ice_var_itd       : convert N-cat to M-cat 
    5353   !!---------------------------------------------------------------------- 
    5454   USE dom_oce        ! ocean space and time domain 
     
    115115      ! 
    116116      ato_i(:,:) = 1._wp - at_i(:,:)         ! open water fraction   
    117       ! 
    118       ALLOCATE( z1_at_i(jpi,jpj) ) 
    119       WHERE( at_i(:,:) > epsi20 )   ;   z1_at_i(:,:) = 1._wp / at_i(:,:) 
    120       ELSEWHERE                     ;   z1_at_i(:,:) = 0._wp 
    121       END WHERE 
    122       ! 
    123       tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:) 
    124       WHERE( at_i(:,:)<=epsi20 ); tm_su(:,:) = rt0; END WHERE       
    125       ! 
     117 
    126118      ! The following fields are calculated for diagnostics and outputs only 
    127119      ! ==> Do not use them for other purposes 
    128120      IF( kn > 1 ) THEN 
    129121         ! 
    130          ALLOCATE( z1_vt_i(jpi,jpj) , z1_vt_s(jpi,jpj) ) 
     122         ALLOCATE( z1_at_i(jpi,jpj) , z1_vt_i(jpi,jpj) , z1_vt_s(jpi,jpj) ) 
     123         WHERE( at_i(:,:) > epsi20 )   ;   z1_at_i(:,:) = 1._wp / at_i(:,:) 
     124         ELSEWHERE                     ;   z1_at_i(:,:) = 0._wp 
     125         END WHERE 
    131126         WHERE( vt_i(:,:) > epsi20 )   ;   z1_vt_i(:,:) = 1._wp / vt_i(:,:) 
    132127         ELSEWHERE                     ;   z1_vt_i(:,:) = 0._wp 
     
    141136         !          
    142137         !                          ! mean temperature (K), salinity and age 
     138         tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:) 
    143139         tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:) 
    144140         om_i (:,:) = SUM( oa_i(:,:,:)              , dim=3 ) * z1_at_i(:,:) 
     
    158154         !                           ! put rt0 where there is no ice 
    159155         WHERE( at_i(:,:)<=epsi20 ) 
     156            tm_su(:,:) = rt0 
    160157            tm_si(:,:) = rt0 
    161158            tm_i (:,:) = rt0 
    162159            tm_s (:,:) = rt0 
    163160         END WHERE 
    164  
    165          DEALLOCATE( z1_vt_i , z1_vt_s ) 
     161         ! 
     162         !                           ! mean melt pond depth 
     163         WHERE( at_ip(:,:) > epsi20 )   ;   hm_ip(:,:) = vt_ip(:,:) / at_ip(:,:) 
     164         ELSEWHERE                      ;   hm_ip(:,:) = 0._wp 
     165         END WHERE          
     166         ! 
     167         DEALLOCATE( z1_at_i , z1_vt_i , z1_vt_s ) 
     168         ! 
    166169      ENDIF 
    167       ! 
    168       DEALLOCATE( z1_at_i ) 
    169170      ! 
    170171   END SUBROUTINE ice_var_agg 
     
    664665      WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0 
    665666      WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0 
    666       IF ( ln_pnd_H12 ) THEN 
     667      IF( ln_pnd_H12 ) THEN 
    667668         WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0 
    668669         WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0 
     
    785786   !! ** Purpose :  converting N-cat ice to jpl ice categories 
    786787   !!------------------------------------------------------------------- 
    787    SUBROUTINE ice_var_itd_1c1c( phti, phts, pati ,       ph_i, ph_s, pa_i, & 
    788       &                         ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i ) 
     788   SUBROUTINE ice_var_itd_1c1c( phti, phts, pati ,                       ph_i, ph_s, pa_i, & 
     789      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 
    789790      !!------------------------------------------------------------------- 
    790791      !! ** Purpose :  converting 1-cat ice to 1 ice category 
     
    792793      REAL(wp), DIMENSION(:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables 
    793794      REAL(wp), DIMENSION(:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables 
    794       REAL(wp), DIMENSION(:), INTENT(in)   , OPTIONAL ::   ptmi, ptms, ptmsu, psmi    ! input  ice/snow temp & sal 
    795       REAL(wp), DIMENSION(:), INTENT(inout), OPTIONAL ::   pt_i, pt_s, pt_su, ps_i    ! output ice/snow temp & sal 
     795      REAL(wp), DIMENSION(:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds 
     796      REAL(wp), DIMENSION(:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds 
    796797      !!------------------------------------------------------------------- 
    797798      ! == thickness and concentration == ! 
     
    800801      pa_i(:) = pati(:) 
    801802      ! 
    802       ! == temperature and salinity == ! 
    803       IF( PRESENT( pt_i  ) )   pt_i (:) = ptmi (:) 
    804       IF( PRESENT( pt_s  ) )   pt_s (:) = ptms (:) 
    805       IF( PRESENT( pt_su ) )   pt_su(:) = ptmsu(:) 
    806       IF( PRESENT( ps_i  ) )   ps_i (:) = psmi (:) 
     803      ! == temperature and salinity and ponds == ! 
     804      pt_i (:) = ptmi (:) 
     805      pt_s (:) = ptms (:) 
     806      pt_su(:) = ptmsu(:) 
     807      ps_i (:) = psmi (:) 
     808      pa_ip(:) = patip(:) 
     809      ph_ip(:) = phtip(:) 
    807810       
    808811   END SUBROUTINE ice_var_itd_1c1c 
    809812 
    810    SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati ,       ph_i, ph_s, pa_i, & 
    811       &                         ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i ) 
     813   SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati ,                       ph_i, ph_s, pa_i, & 
     814      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 
    812815      !!------------------------------------------------------------------- 
    813816      !! ** Purpose :  converting N-cat ice to 1 ice category 
     
    815818      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables 
    816819      REAL(wp), DIMENSION(:)  , INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables 
    817       REAL(wp), DIMENSION(:,:), INTENT(in)   , OPTIONAL ::   ptmi, ptms, ptmsu, psmi    ! input  ice/snow temp & sal 
    818       REAL(wp), DIMENSION(:)  , INTENT(inout), OPTIONAL ::   pt_i, pt_s, pt_su, ps_i    ! output ice/snow temp & sal 
     820      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds 
     821      REAL(wp), DIMENSION(:)  , INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds 
    819822      ! 
    820823      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z1_ai, z1_vi, z1_vs 
     
    826829      ! 
    827830      ! == thickness and concentration == ! 
    828       ALLOCATE( z1_ai(idim) ) 
     831      ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim) ) 
    829832      ! 
    830833      pa_i(:) = SUM( pati(:,:), dim=2 ) 
     
    838841      ! 
    839842      ! == temperature and salinity == ! 
    840       IF( PRESENT( pt_i ) .OR. PRESENT( pt_s ) .OR. PRESENT( pt_su ) .OR. PRESENT( ps_i ) ) THEN 
    841          ! 
    842          ALLOCATE( z1_vi(idim), z1_vs(idim) ) 
    843          ! 
    844          WHERE( ( pa_i(:) * ph_i(:) ) /= 0._wp )   ;   z1_vi(:) = 1._wp / ( pa_i(:) * ph_i(:) ) 
    845          ELSEWHERE                                 ;   z1_vi(:) = 0._wp 
    846          END WHERE 
    847          WHERE( ( pa_i(:) * ph_s(:) ) /= 0._wp )   ;   z1_vs(:) = 1._wp / ( pa_i(:) * ph_s(:) ) 
    848          ELSEWHERE                                 ;   z1_vs(:) = 0._wp 
    849          END WHERE 
    850          ! 
    851          IF( PRESENT( pt_i  ) )   pt_i (:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 
    852          IF( PRESENT( pt_s  ) )   pt_s (:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:) 
    853          IF( PRESENT( pt_su ) )   pt_su(:) = SUM( ptmsu(:,:) * pati(:,:)            , dim=2 ) * z1_ai(:) 
    854          IF( PRESENT( ps_i  ) )   ps_i (:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 
    855          ! 
    856          DEALLOCATE( z1_vi, z1_vs ) 
    857          ! 
    858       ENDIF 
    859       ! 
    860       DEALLOCATE( z1_ai ) 
     843      WHERE( ( pa_i(:) * ph_i(:) ) /= 0._wp )   ;   z1_vi(:) = 1._wp / ( pa_i(:) * ph_i(:) ) 
     844      ELSEWHERE                                 ;   z1_vi(:) = 0._wp 
     845      END WHERE 
     846      WHERE( ( pa_i(:) * ph_s(:) ) /= 0._wp )   ;   z1_vs(:) = 1._wp / ( pa_i(:) * ph_s(:) ) 
     847      ELSEWHERE                                 ;   z1_vs(:) = 0._wp 
     848      END WHERE 
     849      pt_i (:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 
     850      pt_s (:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:) 
     851      pt_su(:) = SUM( ptmsu(:,:) * pati(:,:)            , dim=2 ) * z1_ai(:) 
     852      ps_i (:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 
     853 
     854      ! == ponds == ! 
     855      pa_ip(:) = SUM( patip(:,:), dim=2 ) 
     856      WHERE( pa_ip(:) /= 0._wp )   ;   ph_ip(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / pa_ip(:) 
     857      ELSEWHERE                    ;   ph_ip(:) = 0._wp 
     858      END WHERE 
     859      ! 
     860      DEALLOCATE( z1_ai, z1_vi, z1_vs ) 
    861861      ! 
    862862   END SUBROUTINE ice_var_itd_Nc1c 
    863863    
    864    SUBROUTINE ice_var_itd_1cMc( phti, phts, pati ,       ph_i, ph_s, pa_i, & 
    865       &                         ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i ) 
     864   SUBROUTINE ice_var_itd_1cMc( phti, phts, pati ,                       ph_i, ph_s, pa_i, & 
     865      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 
    866866      !!------------------------------------------------------------------- 
    867867      !! 
     
    894894      REAL(wp), DIMENSION(:),   INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables 
    895895      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables 
    896       REAL(wp), DIMENSION(:)  , INTENT(in)   , OPTIONAL ::   ptmi, ptms, ptmsu, psmi    ! input  ice/snow temp & sal 
    897       REAL(wp), DIMENSION(:,:), INTENT(inout), OPTIONAL ::   pt_i, pt_s, pt_su, ps_i    ! output ice/snow temp & sal 
     896      REAL(wp), DIMENSION(:)  , INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds 
     897      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds 
    898898      ! 
    899899      INTEGER , DIMENSION(4) ::   itest 
     900      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   zfra 
    900901      INTEGER  ::   ji, jk, jl 
    901902      INTEGER  ::   idim, i_fill, jl0   
     
    10101011      ! 
    10111012      ! == temperature and salinity == ! 
    1012       IF( PRESENT( pt_i  ) ) THEN 
    1013          DO jl = 1, jpl 
    1014             pt_i(:,jl) = ptmi (:) 
    1015          END DO 
    1016       ENDIF 
    1017       IF( PRESENT( pt_s  ) ) THEN 
    1018          DO jl = 1, jpl 
    1019             pt_s (:,jl) = ptms (:) 
    1020          END DO 
    1021       ENDIF 
    1022       IF( PRESENT( pt_su ) ) THEN 
    1023          DO jl = 1, jpl 
    1024             pt_su(:,jl) = ptmsu(:) 
    1025          END DO 
    1026       ENDIF 
    1027       IF( PRESENT( ps_i  ) ) THEN 
    1028          DO jl = 1, jpl 
    1029             ps_i (:,jl) = psmi (:) 
    1030          END DO 
    1031       ENDIF 
     1013      DO jl = 1, jpl 
     1014         pt_i (:,jl) = ptmi (:) 
     1015         pt_s (:,jl) = ptms (:) 
     1016         pt_su(:,jl) = ptmsu(:) 
     1017         ps_i (:,jl) = psmi (:) 
     1018         ps_i (:,jl) = psmi (:)          
     1019      END DO 
     1020      ! 
     1021      ! == ponds == ! 
     1022      ALLOCATE( zfra(idim) ) 
     1023      ! keep the same pond fraction atip/ati for each category 
     1024      WHERE( pati(:) /= 0._wp )   ;   zfra(:) = patip(:) / pati(:) 
     1025      ELSEWHERE                   ;   zfra(:) = 0._wp 
     1026      END WHERE 
     1027      DO jl = 1, jpl 
     1028         pa_ip(:,jl) = zfra(:) * pa_i(:,jl) 
     1029      END DO 
     1030      ! keep the same v_ip/v_i ratio for each category 
     1031      WHERE( ( phti(:) * pati(:) ) /= 0._wp )   ;   zfra(:) = ( phtip(:) * patip(:) ) / ( phti(:) * pati(:) ) 
     1032      ELSEWHERE                                 ;   zfra(:) = 0._wp 
     1033      END WHERE 
     1034      DO jl = 1, jpl 
     1035         WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl) 
     1036         ELSEWHERE                       ;   ph_ip(:,jl) = 0._wp 
     1037         END WHERE 
     1038      END DO 
     1039      DEALLOCATE( zfra ) 
    10321040      ! 
    10331041   END SUBROUTINE ice_var_itd_1cMc 
    10341042 
    1035    SUBROUTINE ice_var_itd_NcMc( phti, phts, pati ,       ph_i, ph_s, pa_i, & 
    1036       &                         ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i ) 
     1043   SUBROUTINE ice_var_itd_NcMc( phti, phts, pati ,                       ph_i, ph_s, pa_i, & 
     1044      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 
    10371045      !!------------------------------------------------------------------- 
    10381046      !! 
     
    10651073      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables 
    10661074      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables 
    1067       REAL(wp), DIMENSION(:,:), INTENT(in)   , OPTIONAL ::   ptmi, ptms, ptmsu, psmi    ! input  ice/snow temp & sal 
    1068       REAL(wp), DIMENSION(:,:), INTENT(inout), OPTIONAL ::   pt_i, pt_s, pt_su, ps_i    ! output ice/snow temp & sal 
     1075      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds 
     1076      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds 
    10691077      ! 
    10701078      INTEGER , ALLOCATABLE, DIMENSION(:,:) ::   jlfil, jlfil2 
    10711079      INTEGER , ALLOCATABLE, DIMENSION(:)   ::   jlmax, jlmin 
    1072       REAL(wp), ALLOCATABLE, DIMENSION(:)   ::   z1_ai, z1_vi, z1_vs, ztmp 
     1080      REAL(wp), ALLOCATABLE, DIMENSION(:)   ::   z1_ai, z1_vi, z1_vs, ztmp, zfra 
    10731081      ! 
    10741082      REAL(wp), PARAMETER ::   ztrans = 0.25_wp 
     
    10881096         pa_i(:,:) = pati(:,:) 
    10891097         ! 
    1090          ! == temperature and salinity == ! 
    1091          IF( PRESENT( pt_i  ) )   pt_i (:,:) = ptmi (:,:) 
    1092          IF( PRESENT( pt_s  ) )   pt_s (:,:) = ptms (:,:) 
    1093          IF( PRESENT( pt_su ) )   pt_su(:,:) = ptmsu(:,:) 
    1094          IF( PRESENT( ps_i  ) )   ps_i (:,:) = psmi (:,:) 
     1098         ! == temperature and salinity and ponds == ! 
     1099         pt_i (:,:) = ptmi (:,:) 
     1100         pt_s (:,:) = ptms (:,:) 
     1101         pt_su(:,:) = ptmsu(:,:) 
     1102         ps_i (:,:) = psmi (:,:) 
     1103         pa_ip(:,:) = patip(:,:) 
     1104         ph_ip(:,:) = phtip(:,:) 
    10951105         !                              ! ---------------------- ! 
    10961106      ELSEIF( icat == 1 ) THEN          ! input cat = 1          ! 
    10971107         !                              ! ---------------------- ! 
    1098          CALL  ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1),            ph_i(:,:), ph_s(:,:), pa_i (:,:) ) 
    1099 !!         CALL  ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1),            ph_i(:,:), ph_s(:,:), pa_i (:,:), & 
    1100 !!            &                    ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:) ) 
     1108         CALL  ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1), & 
     1109            &                    ph_i(:,:), ph_s(:,:), pa_i (:,:), & 
     1110            &                    ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), & 
     1111            &                    pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:)  ) 
    11011112         !                              ! ---------------------- ! 
    11021113      ELSEIF( jpl == 1 ) THEN           ! output cat = 1         ! 
    11031114         !                              ! ---------------------- ! 
    1104          CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:),            ph_i(:,1), ph_s(:,1), pa_i (:,1) ) 
    1105 !!         CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:),            ph_i(:,1), ph_s(:,1), pa_i (:,1), & 
    1106 !!            &                    ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1) )          
     1115         CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:), & 
     1116            &                    ph_i(:,1), ph_s(:,1), pa_i (:,1), & 
     1117            &                    ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), & 
     1118            &                    pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1)  ) 
    11071119         !                              ! ----------------------- ! 
    11081120      ELSE                              ! input cat /= output cat ! 
     
    11941206         ! == temperature and salinity == ! 
    11951207         ! 
    1196          IF( PRESENT( pt_i ) .OR. PRESENT( pt_s ) .OR. PRESENT( pt_su ) .OR. PRESENT( ps_i ) ) THEN 
    1197             ! 
    1198             ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim), ztmp(idim) ) 
    1199             ! 
    1200             WHERE( SUM( pa_i(:,:), dim=2 ) /= 0._wp )               ;   z1_ai(:) = 1._wp / SUM( pa_i(:,:), dim=2 ) 
    1201             ELSEWHERE                                               ;   z1_ai(:) = 0._wp 
     1208         ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim), ztmp(idim) ) 
     1209         ! 
     1210         WHERE( SUM( pa_i(:,:), dim=2 ) /= 0._wp )               ;   z1_ai(:) = 1._wp / SUM( pa_i(:,:), dim=2 ) 
     1211         ELSEWHERE                                               ;   z1_ai(:) = 0._wp 
     1212         END WHERE 
     1213         WHERE( SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) /= 0._wp )   ;   z1_vi(:) = 1._wp / SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) 
     1214         ELSEWHERE                                               ;   z1_vi(:) = 0._wp 
     1215         END WHERE 
     1216         WHERE( SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) /= 0._wp )   ;   z1_vs(:) = 1._wp / SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) 
     1217         ELSEWHERE                                               ;   z1_vs(:) = 0._wp 
     1218         END WHERE 
     1219         ! 
     1220         ! fill all the categories with the same value 
     1221         ztmp(:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 
     1222         DO jl = 1, jpl 
     1223            pt_i (:,jl) = ztmp(:) 
     1224         END DO 
     1225         ztmp(:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:) 
     1226         DO jl = 1, jpl 
     1227            pt_s (:,jl) = ztmp(:) 
     1228         END DO 
     1229         ztmp(:) = SUM( ptmsu(:,:) * pati(:,:)            , dim=2 ) * z1_ai(:) 
     1230         DO jl = 1, jpl 
     1231            pt_su(:,jl) = ztmp(:) 
     1232         END DO 
     1233         ztmp(:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 
     1234         DO jl = 1, jpl 
     1235            ps_i (:,jl) = ztmp(:) 
     1236         END DO 
     1237         ! 
     1238         DEALLOCATE( z1_ai, z1_vi, z1_vs, ztmp ) 
     1239         ! 
     1240         ! == ponds == ! 
     1241         ALLOCATE( zfra(idim) ) 
     1242         ! keep the same pond fraction atip/ati for each category 
     1243         WHERE( SUM( pati(:,:), dim=2 ) /= 0._wp )   ;   zfra(:) = SUM( patip(:,:), dim=2 ) / SUM( pati(:,:), dim=2 ) 
     1244         ELSEWHERE                                   ;   zfra(:) = 0._wp 
     1245         END WHERE 
     1246         DO jl = 1, jpl 
     1247            pa_ip(:,jl) = zfra(:) * pa_i(:,jl) 
     1248         END DO 
     1249         ! keep the same v_ip/v_i ratio for each category 
     1250         WHERE( SUM( phti(:,:) * pati(:,:), dim=2 ) /= 0._wp ) 
     1251            zfra(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / SUM( phti(:,:) * pati(:,:), dim=2 ) 
     1252         ELSEWHERE 
     1253            zfra(:) = 0._wp 
     1254         END WHERE 
     1255         DO jl = 1, jpl 
     1256            WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl) 
     1257            ELSEWHERE                       ;   ph_ip(:,jl) = 0._wp 
    12021258            END WHERE 
    1203             WHERE( SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) /= 0._wp )   ;   z1_vi(:) = 1._wp / SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) 
    1204             ELSEWHERE                                               ;   z1_vi(:) = 0._wp 
    1205             END WHERE 
    1206             WHERE( SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) /= 0._wp )   ;   z1_vs(:) = 1._wp / SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) 
    1207             ELSEWHERE                                               ;   z1_vs(:) = 0._wp 
    1208             END WHERE 
    1209             ! 
    1210             ! fill all the categories with the same value 
    1211             IF( PRESENT( pt_i  ) ) THEN 
    1212                ztmp(:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 
    1213                DO jl = 1, jpl 
    1214                   pt_i (:,jl) = ztmp(:) 
    1215                END DO 
    1216             ENDIF 
    1217             IF( PRESENT( pt_s  ) ) THEN 
    1218                ztmp(:) =  SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:) 
    1219                DO jl = 1, jpl 
    1220                   pt_s (:,jl) = ztmp(:) 
    1221                END DO 
    1222             ENDIF 
    1223             IF( PRESENT( pt_su ) ) THEN 
    1224                ztmp(:) = SUM( ptmsu(:,:) * pati(:,:)            , dim=2 ) * z1_ai(:) 
    1225                DO jl = 1, jpl 
    1226                   pt_su(:,jl) = ztmp(:) 
    1227                END DO 
    1228             ENDIF 
    1229             IF( PRESENT( ps_i  ) ) THEN 
    1230                ztmp(:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 
    1231                DO jl = 1, jpl 
    1232                   ps_i (:,jl) = ztmp(:) 
    1233                END DO 
    1234             ENDIF 
    1235             ! 
    1236             DEALLOCATE( z1_ai, z1_vi, z1_vs, ztmp ) 
    1237             ! 
    1238          ENDIF 
     1259         END DO 
     1260         DEALLOCATE( zfra ) 
    12391261         ! 
    12401262      ENDIF 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icewri.F90

    r11413 r11586  
    114114      ! melt ponds 
    115115      IF( iom_use('iceapnd' ) )   CALL iom_put( 'iceapnd', at_ip  * zmsk00      )                                           ! melt pond total fraction 
     116      IF( iom_use('icehpnd' ) )   CALL iom_put( 'icehpnd', hm_ip  * zmsk00      )                                           ! melt pond depth 
    116117      IF( iom_use('icevpnd' ) )   CALL iom_put( 'icevpnd', vt_ip  * zmsk00      )                                           ! melt pond total volume per unit area 
    117118      ! salt 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/BDY/bdy_oce.F90

    r11223 r11586  
    5757      REAL(wp), POINTER, DIMENSION(:,:) ::  h_i    !: Now ice  thickness climatology 
    5858      REAL(wp), POINTER, DIMENSION(:,:) ::  h_s    !: now snow thickness 
     59      REAL(wp), POINTER, DIMENSION(:,:) ::  t_i    !: now ice  temperature 
     60      REAL(wp), POINTER, DIMENSION(:,:) ::  t_s    !: now snow temperature 
     61      REAL(wp), POINTER, DIMENSION(:,:) ::  tsu    !: now surf temperature 
     62      REAL(wp), POINTER, DIMENSION(:,:) ::  s_i    !: now ice  salinity 
     63      REAL(wp), POINTER, DIMENSION(:,:) ::  aip    !: now ice  pond concentration 
     64      REAL(wp), POINTER, DIMENSION(:,:) ::  hip    !: now ice  pond depth 
    5965#if defined key_top 
    6066      CHARACTER(LEN=20)                   :: cn_obc  !: type of boundary condition to apply 
     
    6874   !! Namelist variables 
    6975   !!---------------------------------------------------------------------- 
     76   !                                                   !!** nambdy ** 
    7077   LOGICAL, PUBLIC            ::   ln_bdy                   !: Unstructured Ocean Boundary Condition 
    7178 
     
    101108   INTEGER , DIMENSION(jp_bdy)          ::   nn_ice_dta     !: = 0 use the initial state as bdy dta ;  
    102109                                                            !: = 1 read it in a NetCDF file 
    103    REAL(wp), DIMENSION(jp_bdy) ::   rn_ice_tem              !: choice of the temperature of incoming sea ice 
    104    REAL(wp), DIMENSION(jp_bdy) ::   rn_ice_sal              !: choice of the salinity    of incoming sea ice 
    105    REAL(wp), DIMENSION(jp_bdy) ::   rn_ice_age              !: choice of the age         of incoming sea ice 
     110   !  
     111   !                                                   !!** nambdy_dta ** 
     112   REAL(wp), DIMENSION(jp_bdy) ::   rice_tem                !: temperature of incoming sea ice 
     113   REAL(wp), DIMENSION(jp_bdy) ::   rice_sal                !: salinity    of incoming sea ice 
     114   REAL(wp), DIMENSION(jp_bdy) ::   rice_age                !: age         of incoming sea ice 
     115   REAL(wp), DIMENSION(jp_bdy) ::   rice_apnd               !: pond conc.  of incoming sea ice 
     116   REAL(wp), DIMENSION(jp_bdy) ::   rice_hpnd               !: pond thick. of incoming sea ice 
    106117   ! 
    107     
    108118   !!---------------------------------------------------------------------- 
    109119   !! Global variables 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/BDY/bdydta.F90

    r11413 r11586  
    4343   PUBLIC   bdy_dta_init     ! routine called by nemogcm.F90 
    4444 
    45    INTEGER , PARAMETER ::   jpbdyfld  = 10    ! maximum number of files to read  
     45   INTEGER , PARAMETER ::   jpbdyfld  = 16    ! maximum number of files to read  
    4646   INTEGER , PARAMETER ::   jp_bdyssh = 1     !  
    4747   INTEGER , PARAMETER ::   jp_bdyu2d = 2     !  
     
    5353   INTEGER , PARAMETER ::   jp_bdya_i = 8     !  
    5454   INTEGER , PARAMETER ::   jp_bdyh_i = 9     !  
    55    INTEGER , PARAMETER ::   jp_bdyh_S = 10    !  
     55   INTEGER , PARAMETER ::   jp_bdyh_s = 10    !  
     56   INTEGER , PARAMETER ::   jp_bdyt_i = 11    !  
     57   INTEGER , PARAMETER ::   jp_bdyt_s = 12    !  
     58   INTEGER , PARAMETER ::   jp_bdytsu = 13    !  
     59   INTEGER , PARAMETER ::   jp_bdys_i = 14    !  
     60   INTEGER , PARAMETER ::   jp_bdyaip = 15    !  
     61   INTEGER , PARAMETER ::   jp_bdyhip = 16    !  
    5662#if ! defined key_si3 
    5763   INTEGER , PARAMETER ::   jpl = 1 
    5864#endif 
    59                                                              ! =F => baroclinic velocities in 3D boundary conditions 
     65 
    6066!$AGRIF_DO_NOT_TREAT 
    6167   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:,:), TARGET ::   bf   ! structure of input fields (file informations, fields read) 
     
    181187                        ii = idx_bdy(jbdy)%nbi(ib,igrd) 
    182188                        ij = idx_bdy(jbdy)%nbj(ib,igrd) 
    183                         dta_bdy(jbdy)%a_i (ib,jl) =  a_i(ii,ij,jl) * tmask(ii,ij,1)  
    184                         dta_bdy(jbdy)%h_i (ib,jl) =  h_i(ii,ij,jl) * tmask(ii,ij,1)  
    185                         dta_bdy(jbdy)%h_s (ib,jl) =  h_s(ii,ij,jl) * tmask(ii,ij,1)  
     189                        dta_bdy(jbdy)%a_i(ib,jl) =  a_i (ii,ij,jl) * tmask(ii,ij,1)  
     190                        dta_bdy(jbdy)%h_i(ib,jl) =  h_i (ii,ij,jl) * tmask(ii,ij,1)  
     191                        dta_bdy(jbdy)%h_s(ib,jl) =  h_s (ii,ij,jl) * tmask(ii,ij,1)  
     192                        dta_bdy(jbdy)%t_i(ib,jl) =  SUM(t_i (ii,ij,:,jl)) * r1_nlay_i * tmask(ii,ij,1)  
     193                        dta_bdy(jbdy)%t_s(ib,jl) =  SUM(t_s (ii,ij,:,jl)) * r1_nlay_s * tmask(ii,ij,1) 
     194                        dta_bdy(jbdy)%tsu(ib,jl) =  t_su(ii,ij,jl) * tmask(ii,ij,1)  
     195                        dta_bdy(jbdy)%s_i(ib,jl) =  s_i (ii,ij,jl) * tmask(ii,ij,1) 
     196                        ! melt ponds 
     197                        dta_bdy(jbdy)%aip(ib,jl) =  a_ip(ii,ij,jl) * tmask(ii,ij,1)  
     198                        dta_bdy(jbdy)%hip(ib,jl) =  h_ip(ii,ij,jl) * tmask(ii,ij,1)  
    186199                     END DO 
    187200                  END DO 
     
    280293 
    281294#if defined key_si3 
    282          ! ice: convert N-cat fields (input) into jpl-cat (output) 
    283295         IF( dta_alias%lneed_ice ) THEN 
    284             ipl = SIZE(bf_alias(jp_bdya_i)%fnow, 3) 
     296            ! fill temperature and salinity arrays 
     297            IF( TRIM(bf_alias(jp_bdyt_i)%clrootname) == 'NOT USED' )   bf_alias(jp_bdyt_i)%fnow(:,1,:) = rice_tem (jbdy) 
     298            IF( TRIM(bf_alias(jp_bdyt_s)%clrootname) == 'NOT USED' )   bf_alias(jp_bdyt_s)%fnow(:,1,:) = rice_tem (jbdy) 
     299            IF( TRIM(bf_alias(jp_bdytsu)%clrootname) == 'NOT USED' )   bf_alias(jp_bdytsu)%fnow(:,1,:) = rice_tem (jbdy) 
     300            IF( TRIM(bf_alias(jp_bdys_i)%clrootname) == 'NOT USED' )   bf_alias(jp_bdys_i)%fnow(:,1,:) = rice_sal (jbdy) 
     301            IF( TRIM(bf_alias(jp_bdyaip)%clrootname) == 'NOT USED' )   bf_alias(jp_bdyaip)%fnow(:,1,:) = rice_apnd(jbdy) * & ! rice_apnd is the pond fraction 
     302               &                                                                         bf_alias(jp_bdya_i)%fnow(:,1,:)     !   ( a_ip = rice_apnd * a_i ) 
     303            IF( TRIM(bf_alias(jp_bdyhip)%clrootname) == 'NOT USED' )   bf_alias(jp_bdyhip)%fnow(:,1,:) = rice_hpnd(jbdy) 
     304            ! if T_su is read and not T_i, set T_i = (T_su + T_freeze)/2 
     305            IF( TRIM(bf_alias(jp_bdytsu)%clrootname) /= 'NOT USED' .AND. TRIM(bf_alias(jp_bdyt_i)%clrootname) == 'NOT USED' ) & 
     306               &   bf_alias(jp_bdyt_i)%fnow(:,1,:) = 0.5_wp * ( bf_alias(jp_bdytsu)%fnow(:,1,:) + 271.15 ) 
     307            ! if T_su is read and not T_s, set T_s = T_su 
     308            IF( TRIM(bf_alias(jp_bdytsu)%clrootname) /= 'NOT USED' .AND. TRIM(bf_alias(jp_bdyt_s)%clrootname) == 'NOT USED' ) & 
     309               &   bf_alias(jp_bdyt_s)%fnow(:,1,:) = bf_alias(jp_bdytsu)%fnow(:,1,:) 
     310            ! if T_s is read and not T_su, set T_su = T_s 
     311            IF( TRIM(bf_alias(jp_bdyt_s)%clrootname) /= 'NOT USED' .AND. TRIM(bf_alias(jp_bdytsu)%clrootname) == 'NOT USED' ) & 
     312               &   bf_alias(jp_bdytsu)%fnow(:,1,:) = bf_alias(jp_bdyt_s)%fnow(:,1,:) 
     313            ! if T_s is read and not T_i, set T_i = (T_s + T_freeze)/2 
     314            IF( TRIM(bf_alias(jp_bdyt_s)%clrootname) /= 'NOT USED' .AND. TRIM(bf_alias(jp_bdyt_i)%clrootname) == 'NOT USED' ) & 
     315               &   bf_alias(jp_bdyt_i)%fnow(:,1,:) = 0.5_wp * ( bf_alias(jp_bdyt_s)%fnow(:,1,:) + 271.15 ) 
     316 
     317            ! make sure ponds = 0 if no ponds scheme 
     318            IF ( .NOT.ln_pnd ) THEN 
     319               bf_alias(jp_bdyaip)%fnow(:,1,:) = 0._wp 
     320               bf_alias(jp_bdyhip)%fnow(:,1,:) = 0._wp 
     321            ENDIF 
     322             
     323            ! convert N-cat fields (input) into jpl-cat (output) 
     324            ipl = SIZE(bf_alias(jp_bdya_i)%fnow, 3)             
    285325            IF( ipl /= jpl ) THEN      ! ice: convert N-cat fields (input) into jpl-cat (output) 
    286                CALL ice_var_itd(bf_alias(jp_bdyh_i)%fnow(:,1,:), bf_alias(jp_bdyh_s)%fnow(:,1,:), bf_alias(jp_bdya_i)%fnow(:,1,:), & 
    287                   &              dta_alias%h_i               , dta_alias%h_s               , dta_alias%a_i                 ) 
     326               CALL ice_var_itd( bf_alias(jp_bdyh_i)%fnow(:,1,:), bf_alias(jp_bdyh_s)%fnow(:,1,:), bf_alias(jp_bdya_i)%fnow(:,1,:), & 
     327                  &              dta_alias%h_i                  , dta_alias%h_s                  , dta_alias%a_i                  , & 
     328                  &              bf_alias(jp_bdyt_i)%fnow(:,1,:), bf_alias(jp_bdyt_s)%fnow(:,1,:), & 
     329                  &              bf_alias(jp_bdytsu)%fnow(:,1,:), bf_alias(jp_bdys_i)%fnow(:,1,:), & 
     330                  &              bf_alias(jp_bdyaip)%fnow(:,1,:), bf_alias(jp_bdyhip)%fnow(:,1,:), & 
     331                  &              dta_alias%t_i                  , dta_alias%t_s                  , & 
     332                  &              dta_alias%tsu                  , dta_alias%s_i                  , & 
     333                  &              dta_alias%aip                  , dta_alias%hip ) 
    288334            ENDIF 
    289335         ENDIF 
     
    332378      !                                                         ! =F => baroclinic velocities in 3D boundary data 
    333379      LOGICAL                                ::   ln_zinterp    ! =T => requires a vertical interpolation of the bdydta 
     380      REAL(wp)                               ::   rn_ice_tem, rn_ice_sal, rn_ice_age, rn_ice_apnd, rn_ice_hpnd  
    334381      INTEGER                                ::   ipk,ipl       ! 
    335382      INTEGER                                ::   idvar         ! variable ID 
     
    341388      LOGICAL                                ::   llneed        ! 
    342389      LOGICAL                                ::   llread        ! 
    343       TYPE(FLD_N), DIMENSION(1), TARGET ::   bn_tem, bn_sal, bn_u3d, bn_v3d   ! must be an array to be used with fld_fill 
    344       TYPE(FLD_N), DIMENSION(1), TARGET ::   bn_ssh, bn_u2d, bn_v2d           ! informations about the fields to be read 
    345       TYPE(FLD_N), DIMENSION(1), TARGET ::   bn_a_i, bn_h_i, bn_h_s       
     390      TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_tem, bn_sal, bn_u3d, bn_v3d   ! must be an array to be used with fld_fill 
     391      TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_ssh, bn_u2d, bn_v2d           ! informations about the fields to be read 
     392      TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_a_i, bn_h_i, bn_h_s, bn_t_i, bn_t_s, bn_tsu, bn_s_i, bn_aip, bn_hip        
    346393      TYPE(FLD_N), DIMENSION(:), POINTER ::   bn_alias                        ! must be an array to be used with fld_fill 
    347394      TYPE(FLD  ), DIMENSION(:), POINTER ::   bf_alias 
    348395      ! 
    349396      NAMELIST/nambdy_dta/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d  
    350       NAMELIST/nambdy_dta/ bn_a_i, bn_h_i, bn_h_s 
     397      NAMELIST/nambdy_dta/ bn_a_i, bn_h_i, bn_h_s, bn_t_i, bn_t_s, bn_tsu, bn_s_i, bn_aip, bn_hip 
     398      NAMELIST/nambdy_dta/ rn_ice_tem, rn_ice_sal, rn_ice_age, rn_ice_apnd, rn_ice_hpnd 
    351399      NAMELIST/nambdy_dta/ ln_full_vel, ln_zinterp 
    352400      !!--------------------------------------------------------------------------- 
     
    402450         ENDIF 
    403451 
     452#if defined key_si3 
     453         IF( .NOT.ln_pnd ) THEN 
     454            rn_ice_apnd = 0. ; rn_ice_hpnd = 0. 
     455            CALL ctl_warn( 'rn_ice_apnd & rn_ice_hpnd = 0 when no ponds' ) 
     456         ENDIF 
     457#endif 
     458 
     459         ! temp, salt, age and ponds of incoming ice 
     460         rice_tem (jbdy) = rn_ice_tem 
     461         rice_sal (jbdy) = rn_ice_sal 
     462         rice_age (jbdy) = rn_ice_age 
     463         rice_apnd(jbdy) = rn_ice_apnd 
     464         rice_hpnd(jbdy) = rn_ice_hpnd 
     465          
     466          
    404467         DO jfld = 1, jpbdyfld 
    405468 
     
    497560            !          ice 
    498561            ! ===================== 
     562            IF(  jfld == jp_bdya_i .OR. jfld == jp_bdyh_i .OR. jfld == jp_bdyh_s .OR. & 
     563               & jfld == jp_bdyt_i .OR. jfld == jp_bdyt_s .OR. jfld == jp_bdytsu .OR. & 
     564               & jfld == jp_bdys_i .OR. jfld == jp_bdyaip .OR. jfld == jp_bdyhip      ) THEN 
     565               igrd = 1                                                    ! T point 
     566               ipk = ipl                                                   ! jpl-cat data 
     567               llneed = dta_bdy(jbdy)%lneed_ice                            ! ice will be needed 
     568               llread = nn_ice_dta(jbdy) == 1                              ! get data from NetCDF file 
     569               iszdim = idx_bdy(jbdy)%nblen(igrd)                          ! length of this bdy on this MPI processus 
     570            ENDIF 
    499571            IF( jfld == jp_bdya_i ) THEN 
    500572               cl3 = 'a_i' 
    501                igrd = 1                                                    ! T point 
    502                ipk = ipl                                                   !  
    503                llneed = dta_bdy(jbdy)%lneed_ice                            ! dta_bdy(jbdy)%a_i will be needed 
    504                llread = nn_ice_dta(jbdy) == 1                              ! get data from NetCDF file 
    505                bf_alias => bf(jp_bdya_i,jbdy:jbdy)                         ! alias for ssh structure of bdy number jbdy 
    506                bn_alias => bn_a_i                                          ! alias for ssh structure of nambdy_dta  
    507                iszdim = idx_bdy(jbdy)%nblen(igrd)                          ! length of this bdy on this MPI processus 
    508            ENDIF 
     573               bf_alias => bf(jp_bdya_i,jbdy:jbdy)                         ! alias for a_i structure of bdy number jbdy 
     574               bn_alias => bn_a_i                                          ! alias for a_i structure of nambdy_dta  
     575            ENDIF 
    509576            IF( jfld == jp_bdyh_i ) THEN 
    510577               cl3 = 'h_i' 
    511                igrd = 1                                                    ! T point 
    512                ipk = ipl                                                   !  
    513                llneed = dta_bdy(jbdy)%lneed_ice                            ! dta_bdy(jbdy)%h_i will be needed 
    514                llread = nn_ice_dta(jbdy) == 1                              ! get data from NetCDF file 
    515                bf_alias => bf(jp_bdyh_i,jbdy:jbdy)                         ! alias for ssh structure of bdy number jbdy 
    516                bn_alias => bn_h_i                                          ! alias for ssh structure of nambdy_dta  
    517                iszdim = idx_bdy(jbdy)%nblen(igrd)                          ! length of this bdy on this MPI processus 
     578               bf_alias => bf(jp_bdyh_i,jbdy:jbdy)                         ! alias for h_i structure of bdy number jbdy 
     579               bn_alias => bn_h_i                                          ! alias for h_i structure of nambdy_dta  
    518580            ENDIF 
    519581            IF( jfld == jp_bdyh_s ) THEN 
    520582               cl3 = 'h_s' 
    521                igrd = 1                                                    ! T point 
    522                ipk = ipl                                                   !  
    523                llneed = dta_bdy(jbdy)%lneed_ice                            ! dta_bdy(jbdy)%h_s will be needed 
    524                llread = nn_ice_dta(jbdy) == 1                              ! get data from NetCDF file 
    525                bf_alias => bf(jp_bdyh_s,jbdy:jbdy)                         ! alias for ssh structure of bdy number jbdy 
    526                bn_alias => bn_h_s                                          ! alias for ssh structure of nambdy_dta  
    527                iszdim = idx_bdy(jbdy)%nblen(igrd)                          ! length of this bdy on this MPI processus 
     583               bf_alias => bf(jp_bdyh_s,jbdy:jbdy)                         ! alias for h_s structure of bdy number jbdy 
     584               bn_alias => bn_h_s                                          ! alias for h_s structure of nambdy_dta  
     585            ENDIF 
     586            IF( jfld == jp_bdyt_i ) THEN 
     587               cl3 = 't_i' 
     588               bf_alias => bf(jp_bdyt_i,jbdy:jbdy)                         ! alias for t_i structure of bdy number jbdy 
     589               bn_alias => bn_t_i                                          ! alias for t_i structure of nambdy_dta  
     590            ENDIF 
     591            IF( jfld == jp_bdyt_s ) THEN 
     592               cl3 = 't_s' 
     593               bf_alias => bf(jp_bdyt_s,jbdy:jbdy)                         ! alias for t_s structure of bdy number jbdy 
     594               bn_alias => bn_t_s                                          ! alias for t_s structure of nambdy_dta  
     595            ENDIF 
     596            IF( jfld == jp_bdytsu ) THEN 
     597               cl3 = 'tsu' 
     598               bf_alias => bf(jp_bdytsu,jbdy:jbdy)                         ! alias for tsu structure of bdy number jbdy 
     599               bn_alias => bn_tsu                                          ! alias for tsu structure of nambdy_dta  
     600            ENDIF 
     601            IF( jfld == jp_bdys_i ) THEN 
     602               cl3 = 's_i' 
     603               bf_alias => bf(jp_bdys_i,jbdy:jbdy)                         ! alias for s_i structure of bdy number jbdy 
     604               bn_alias => bn_s_i                                          ! alias for s_i structure of nambdy_dta  
     605            ENDIF 
     606            IF( jfld == jp_bdyaip ) THEN 
     607               cl3 = 'aip' 
     608               bf_alias => bf(jp_bdyaip,jbdy:jbdy)                         ! alias for aip structure of bdy number jbdy 
     609               bn_alias => bn_aip                                          ! alias for aip structure of nambdy_dta  
     610            ENDIF 
     611            IF( jfld == jp_bdyhip ) THEN 
     612               cl3 = 'hip' 
     613               bf_alias => bf(jp_bdyhip,jbdy:jbdy)                         ! alias for hip structure of bdy number jbdy 
     614               bn_alias => bn_hip                                          ! alias for hip structure of nambdy_dta  
    528615            ENDIF 
    529616 
     
    542629 
    543630               ! associate the pointer and get rid of the dimensions with a size equal to 1 
    544                IF( jfld == jp_bdyssh           ) dta_bdy(jbdy)%ssh => bf_alias(1)%fnow(:,1,1) 
    545                IF( jfld == jp_bdyu2d           ) dta_bdy(jbdy)%u2d => bf_alias(1)%fnow(:,1,1) 
    546                IF( jfld == jp_bdyv2d           ) dta_bdy(jbdy)%v2d => bf_alias(1)%fnow(:,1,1) 
    547                IF( jfld == jp_bdyu3d           ) dta_bdy(jbdy)%u3d => bf_alias(1)%fnow(:,1,:) 
    548                IF( jfld == jp_bdyv3d           ) dta_bdy(jbdy)%v3d => bf_alias(1)%fnow(:,1,:) 
    549                IF( jfld == jp_bdytem           ) dta_bdy(jbdy)%tem => bf_alias(1)%fnow(:,1,:) 
    550                IF( jfld == jp_bdysal           ) dta_bdy(jbdy)%sal => bf_alias(1)%fnow(:,1,:) 
     631               IF( jfld == jp_bdyssh )        dta_bdy(jbdy)%ssh => bf_alias(1)%fnow(:,1,1) 
     632               IF( jfld == jp_bdyu2d )        dta_bdy(jbdy)%u2d => bf_alias(1)%fnow(:,1,1) 
     633               IF( jfld == jp_bdyv2d )        dta_bdy(jbdy)%v2d => bf_alias(1)%fnow(:,1,1) 
     634               IF( jfld == jp_bdyu3d )        dta_bdy(jbdy)%u3d => bf_alias(1)%fnow(:,1,:) 
     635               IF( jfld == jp_bdyv3d )        dta_bdy(jbdy)%v3d => bf_alias(1)%fnow(:,1,:) 
     636               IF( jfld == jp_bdytem )        dta_bdy(jbdy)%tem => bf_alias(1)%fnow(:,1,:) 
     637               IF( jfld == jp_bdysal )        dta_bdy(jbdy)%sal => bf_alias(1)%fnow(:,1,:) 
    551638               IF( jfld == jp_bdya_i ) THEN 
    552639                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%a_i => bf_alias(1)%fnow(:,1,:) 
     
    562649                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%h_s => bf_alias(1)%fnow(:,1,:) 
    563650                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%h_s(iszdim,jpl) ) 
     651                  ENDIF 
     652               ENDIF 
     653               IF( jfld == jp_bdyt_i ) THEN 
     654                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%t_i => bf_alias(1)%fnow(:,1,:) 
     655                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%t_i(iszdim,jpl) ) 
     656                  ENDIF 
     657               ENDIF 
     658               IF( jfld == jp_bdyt_s ) THEN 
     659                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%t_s => bf_alias(1)%fnow(:,1,:) 
     660                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%t_s(iszdim,jpl) ) 
     661                  ENDIF 
     662               ENDIF 
     663               IF( jfld == jp_bdytsu ) THEN 
     664                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%tsu => bf_alias(1)%fnow(:,1,:) 
     665                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%tsu(iszdim,jpl) ) 
     666                  ENDIF 
     667               ENDIF 
     668               IF( jfld == jp_bdys_i ) THEN 
     669                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%s_i => bf_alias(1)%fnow(:,1,:) 
     670                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%s_i(iszdim,jpl) ) 
     671                  ENDIF 
     672               ENDIF 
     673               IF( jfld == jp_bdyaip ) THEN 
     674                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%aip => bf_alias(1)%fnow(:,1,:) 
     675                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%aip(iszdim,jpl) ) 
     676                  ENDIF 
     677               ENDIF 
     678               IF( jfld == jp_bdyhip ) THEN 
     679                  IF( ipk == jpl ) THEN   ;   dta_bdy(jbdy)%hip => bf_alias(1)%fnow(:,1,:) 
     680                  ELSE                    ;   ALLOCATE( dta_bdy(jbdy)%hip(iszdim,jpl) ) 
    564681                  ENDIF 
    565682               ENDIF 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/BDY/bdyice.F90

    r11210 r11586  
    6363      IF( ln_timing    )   CALL timing_start('bdy_ice_thd')                                                            ! timing 
    6464      IF( ln_icediachk )   CALL ice_cons_hsm(0,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
     65      IF( ln_icediachk )   CALL ice_cons2D  (0,'bdy_ice_thd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation 
    6566      ! 
    6667      CALL ice_var_glo2eqv 
     
    109110      ! 
    110111      ! controls 
     112      IF( ln_icectl    )   CALL ice_prt     ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )                        ! prints 
    111113      IF( ln_icediachk )   CALL ice_cons_hsm(1,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    112       IF( ln_icectl    )   CALL ice_prt     ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )                        ! prints 
     114      IF( ln_icediachk )   CALL ice_cons2D  (1,'bdy_ice_thd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation 
    113115      IF( ln_timing    )   CALL timing_stop ('bdy_ice_thd')                                                            ! timing 
    114116      ! 
     
    152154            zwgt  = idx%nbw(i_bdy,jgrd) 
    153155            zwgt1 = 1.e0 - idx%nbw(i_bdy,jgrd) 
    154             a_i(ji,jj,jl) = ( a_i(ji,jj,jl) * zwgt1 + dta%a_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Leads fraction  
    155             h_i(ji,jj,jl) = ( h_i(ji,jj,jl) * zwgt1 + dta%h_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice depth  
    156             h_s(ji,jj,jl) = ( h_s(ji,jj,jl) * zwgt1 + dta%h_s(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Snow depth 
    157  
     156            a_i (ji,jj,  jl) = ( a_i (ji,jj,  jl) * zwgt1 + dta%a_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  concentration  
     157            h_i (ji,jj,  jl) = ( h_i (ji,jj,  jl) * zwgt1 + dta%h_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  depth  
     158            h_s (ji,jj,  jl) = ( h_s (ji,jj,  jl) * zwgt1 + dta%h_s(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Snow depth 
     159            t_i (ji,jj,:,jl) = ( t_i (ji,jj,:,jl) * zwgt1 + dta%t_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  temperature 
     160            t_s (ji,jj,:,jl) = ( t_s (ji,jj,:,jl) * zwgt1 + dta%t_s(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Snow temperature 
     161            t_su(ji,jj,  jl) = ( t_su(ji,jj,  jl) * zwgt1 + dta%tsu(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Surf temperature 
     162            s_i (ji,jj,  jl) = ( s_i (ji,jj,  jl) * zwgt1 + dta%s_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  salinity 
     163            a_ip(ji,jj,  jl) = ( a_ip(ji,jj,  jl) * zwgt1 + dta%aip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  pond concentration 
     164            h_ip(ji,jj,  jl) = ( h_ip(ji,jj,  jl) * zwgt1 + dta%hip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  pond depth 
     165            ! 
     166            sz_i(ji,jj,:,jl) = s_i(ji,jj,jl) 
     167            ! 
     168            ! make sure ponds = 0 if no ponds scheme 
     169            IF( .NOT.ln_pnd ) THEN 
     170               a_ip(ji,jj,jl) = 0._wp 
     171               h_ip(ji,jj,jl) = 0._wp 
     172            ENDIF 
     173            ! 
    158174            ! ----------------- 
    159175            ! Pathological case 
     
    170186            h_i(ji,jj,jl) = MIN( hi_max(jl), h_i(ji,jj,jl) + zdh ) 
    171187            h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoi / rhos )  
    172  
     188            ! 
    173189         ENDDO 
    174190      ENDDO 
     
    206222            IF( a_i(ib,jb,jl) > 0._wp ) THEN   ! there is ice at the boundary 
    207223               ! 
    208                a_i(ji,jj,jl) = a_i(ib,jb,jl) ! concentration 
    209                h_i(ji,jj,jl) = h_i(ib,jb,jl) ! thickness ice 
    210                h_s(ji,jj,jl) = h_s(ib,jb,jl) ! thickness snw 
    211                ! 
    212                SELECT CASE( jpbound ) 
    213                   ! 
    214                CASE( 0 )   ! velocity is inward 
    215                   ! 
    216                   oa_i(ji,jj,  jl) = rn_ice_age(jbdy) * a_i(ji,jj,jl) ! age 
    217                   a_ip(ji,jj,  jl) = 0._wp                            ! pond concentration 
    218                   v_ip(ji,jj,  jl) = 0._wp                            ! pond volume 
    219                   t_su(ji,jj,  jl) = rn_ice_tem(jbdy)                 ! temperature surface 
    220                   t_s (ji,jj,:,jl) = rn_ice_tem(jbdy)                 ! temperature snw 
    221                   t_i (ji,jj,:,jl) = rn_ice_tem(jbdy)                 ! temperature ice 
    222                   s_i (ji,jj,  jl) = rn_ice_sal(jbdy)                 ! salinity 
    223                   sz_i(ji,jj,:,jl) = rn_ice_sal(jbdy)                 ! salinity profile 
    224                   ! 
    225                CASE( 1 )   ! velocity is outward 
    226                   ! 
    227                   oa_i(ji,jj,  jl) = oa_i(ib,jb,  jl) ! age 
    228                   a_ip(ji,jj,  jl) = a_ip(ib,jb,  jl) ! pond concentration 
    229                   v_ip(ji,jj,  jl) = v_ip(ib,jb,  jl) ! pond volume 
    230                   t_su(ji,jj,  jl) = t_su(ib,jb,  jl) ! temperature surface 
    231                   t_s (ji,jj,:,jl) = t_s (ib,jb,:,jl) ! temperature snw 
    232                   t_i (ji,jj,:,jl) = t_i (ib,jb,:,jl) ! temperature ice 
    233                   s_i (ji,jj,  jl) = s_i (ib,jb,  jl) ! salinity 
    234                   sz_i(ji,jj,:,jl) = sz_i(ib,jb,:,jl) ! salinity profile 
    235                   ! 
    236                END SELECT 
     224               a_i (ji,jj,  jl) = a_i (ib,jb,  jl) 
     225               h_i (ji,jj,  jl) = h_i (ib,jb,  jl) 
     226               h_s (ji,jj,  jl) = h_s (ib,jb,  jl) 
     227               t_i (ji,jj,:,jl) = t_i (ib,jb,:,jl) 
     228               t_s (ji,jj,:,jl) = t_s (ib,jb,:,jl) 
     229               t_su(ji,jj,  jl) = t_su(ib,jb,  jl) 
     230               s_i (ji,jj,  jl) = s_i (ib,jb,  jl) 
     231               a_ip(ji,jj,  jl) = a_ip(ib,jb,  jl) 
     232               h_ip(ji,jj,  jl) = h_ip(ib,jb,  jl) 
     233               ! 
     234               sz_i(ji,jj,:,jl) = sz_i(ib,jb,:,jl) 
     235               ! 
     236               ! ice age 
     237               IF    ( jpbound == 0 ) THEN  ! velocity is inward 
     238                  oa_i(ji,jj,jl) = rice_age(jbdy) * a_i(ji,jj,jl) 
     239               ELSEIF( jpbound == 1 ) THEN  ! velocity is outward 
     240                  oa_i(ji,jj,jl) = oa_i(ib,jb,jl) 
     241               ENDIF 
    237242               ! 
    238243               IF( nn_icesal == 1 ) THEN     ! if constant salinity 
     
    259264               END DO 
    260265               ! 
     266               ! melt ponds 
     267               IF( a_i(ji,jj,jl) > epsi10 ) THEN 
     268                  a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i (ji,jj,jl) 
     269               ELSE 
     270                  a_ip_frac(ji,jj,jl) = 0._wp 
     271               ENDIF 
     272               v_ip(ji,jj,jl) = h_ip(ji,jj,jl) * a_ip(ji,jj,jl) 
     273               ! 
    261274            ELSE   ! no ice at the boundary 
    262275               ! 
     
    270283               t_s (ji,jj,:,jl) = rt0 
    271284               t_i (ji,jj,:,jl) = rt0  
     285 
     286               a_ip_frac(ji,jj,jl) = 0._wp 
     287               h_ip     (ji,jj,jl) = 0._wp 
     288               a_ip     (ji,jj,jl) = 0._wp 
     289               v_ip     (ji,jj,jl) = 0._wp 
    272290                
    273291               IF( nn_icesal == 1 ) THEN     ! if constant salinity 
     
    372390                     jj    = idx_bdy(jbdy)%nbj(i_bdy,jgrd) 
    373391                     zflag = idx_bdy(jbdy)%flagv(i_bdy,jgrd) 
    374                      !    ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨     !  ¨¨¨¨ïce¨¨¨(jj+1)¨¨     ! ¨¨¨¨¨¨ö¨¨¨¨(jj+1)        
     392                     !                         !      ice   (jj+1)       !       o    (jj+1) 
    375393                     !       ^    (jj  )       !       ^    (jj  )       !       ^    (jj  )        
    376394                     !      ice   (jj  )       !       o    (jj  )       !       o    (jj  )        
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/BDY/bdyini.F90

    r11413 r11586  
    6767         &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, & 
    6868         &             cn_ice, nn_ice_dta,                                     & 
    69          &             rn_ice_tem, rn_ice_sal, rn_ice_age,                     & 
    7069         &             ln_vol, nn_volctl, nn_rimwidth 
    7170         ! 
     
    9493      cn_ice         (2:jp_bdy) = cn_ice         (1) 
    9594      nn_ice_dta     (2:jp_bdy) = nn_ice_dta     (1) 
    96       rn_ice_tem     (2:jp_bdy) = rn_ice_tem     (1) 
    97       rn_ice_sal     (2:jp_bdy) = rn_ice_sal     (1) 
    98       rn_ice_age     (2:jp_bdy) = rn_ice_age     (1) 
    9995      REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist :Unstructured open boundaries 
    10096      READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 902 ) 
     
    328324            CASE DEFAULT   ;   CALL ctl_stop( 'nn_ice_dta must be 0 or 1' ) 
    329325            END SELECT 
    330             WRITE(numout,*) 
    331             WRITE(numout,*) '      tem of bdy sea-ice = ', rn_ice_tem(ib_bdy)          
    332             WRITE(numout,*) '      sal of bdy sea-ice = ', rn_ice_sal(ib_bdy)          
    333             WRITE(numout,*) '      age of bdy sea-ice = ', rn_ice_age(ib_bdy)          
    334326         ENDIF 
    335327#else 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DIA/diacfl.F90

    r10824 r11586  
    2929   REAL(wp)              ::   rCu_max, rCv_max, rCw_max   ! associated run max Courant number  
    3030 
    31 !!gm CAUTION: need to declare these arrays here, otherwise the calculation fails in multi-proc ! 
    32 !!gm          I don't understand why. 
    33    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zCu_cfl, zCv_cfl, zCw_cfl         ! workspace 
    34 !!gm end 
    35  
    3631   PUBLIC   dia_cfl       ! routine called by step.F90 
    3732   PUBLIC   dia_cfl_init  ! routine called by nemogcm 
     
    5550      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    5651      ! 
    57       INTEGER                ::   ji, jj, jk                            ! dummy loop indices 
    58       REAL(wp)               ::   z2dt, zCu_max, zCv_max, zCw_max       ! local scalars 
    59       INTEGER , DIMENSION(3) ::   iloc_u , iloc_v , iloc_w , iloc       ! workspace 
    60 !!gm this does not work      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zCu_cfl, zCv_cfl, zCw_cfl         ! workspace 
     52      INTEGER                          ::   ji, jj, jk                       ! dummy loop indices 
     53      REAL(wp)                         ::   z2dt, zCu_max, zCv_max, zCw_max  ! local scalars 
     54      INTEGER , DIMENSION(3)           ::   iloc_u , iloc_v , iloc_w , iloc  ! workspace 
     55      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zCu_cfl, zCv_cfl, zCw_cfl        ! workspace 
    6156      !!---------------------------------------------------------------------- 
    6257      ! 
     
    7166      DO jk = 1, jpk       ! calculate Courant numbers 
    7267         DO jj = 1, jpj 
    73             DO ji = 1, fs_jpim1   ! vector opt. 
     68            DO ji = 1, jpi 
    7469               zCu_cfl(ji,jj,jk) = ABS( un(ji,jj,jk) ) * z2dt / e1u  (ji,jj)      ! for i-direction 
    7570               zCv_cfl(ji,jj,jk) = ABS( vn(ji,jj,jk) ) * z2dt / e2v  (ji,jj)      ! for j-direction 
     
    111106      !                    ! write out to file 
    112107      IF( lwp ) THEN 
    113          WRITE(numcfl,FMT='(2x,i4,5x,a6,4x,f7.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zCu_max, iloc_u(1), iloc_u(2), iloc_u(3) 
     108         WRITE(numcfl,FMT='(2x,i6,3x,a6,4x,f7.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zCu_max, iloc_u(1), iloc_u(2), iloc_u(3) 
    114109         WRITE(numcfl,FMT='(11x,     a6,4x,f7.4,1x,i4,1x,i4,1x,i4)')     'Max Cv', zCv_max, iloc_v(1), iloc_v(2), iloc_v(3) 
    115110         WRITE(numcfl,FMT='(11x,     a6,4x,f7.4,1x,i4,1x,i4,1x,i4)')     'Max Cw', zCw_max, iloc_w(1), iloc_w(2), iloc_w(3) 
     
    172167      rCw_max = 0._wp 
    173168      ! 
    174 !!gm required to work 
    175       ALLOCATE ( zCu_cfl(jpi,jpj,jpk), zCv_cfl(jpi,jpj,jpk), zCw_cfl(jpi,jpj,jpk) ) 
    176 !!gm end 
    177       !       
    178169   END SUBROUTINE dia_cfl_init 
    179170 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DIA/diawri.F90

    r11413 r11586  
    7171#endif 
    7272   INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file 
    73    INTEGER ::          nb_T, ndim_bT                 ! grid_T file 
     73   INTEGER ::          nb_T              , ndim_bT   ! grid_T file 
    7474   INTEGER ::   nid_U, nz_U, nh_U, ndim_U, ndim_hU   ! grid_U file 
    7575   INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file 
    7676   INTEGER ::   nid_W, nz_W, nh_W                    ! grid_W file 
    77    INTEGER ::   ndim_A, ndim_hA                      ! ABL file    
    78    INTEGER ::   nid_A, nz_A, nh_A                    ! grid_ABL file    
     77   INTEGER ::   nid_A, nz_A, nh_A, ndim_A, ndim_hA   ! grid_ABL file    
    7978   INTEGER ::   ndex(1)                              ! ??? 
    8079   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV 
     
    216215      ENDIF 
    217216 
     217      IF( ln_zad_Aimp ) wn = wn + wi               ! Recombine explicit and implicit parts of vertical velocity for diagnostic output 
     218      ! 
    218219      CALL iom_put( "woce", wn )                   ! vertical velocity 
    219220      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value 
     
    226227         IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) 
    227228      ENDIF 
     229      ! 
     230      IF( ln_zad_Aimp ) wn = wn - wi               ! Remove implicit part of vertical velocity that was added for diagnostic output 
    228231 
    229232      CALL iom_put( "avt" , avt )                  ! T vert. eddy diff. coef. 
     
    415418         &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     & 
    416419         &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) ) 
     420         ! 
    417421     dia_wri_alloc = MAXVAL(ierr) 
    418422      CALL mpp_sum( 'diawri', dia_wri_alloc ) 
     
    460464         ninist = 0 
    461465      ENDIF 
    462       
    463466      ! 
    464467      IF( nn_write == -1 )   RETURN   ! we will never do any output 
     
    488491      ijmi = 1      ;      ijma = jpj 
    489492      ipk = jpk 
    490      IF(ln_abl) ipka = jpkam1 
     493      IF(ln_abl) ipka = jpkam1 
    491494 
    492495      ! define time axis 
     
    724727        
    725728         clmx ="l_max(only(x))"    ! max index on a period 
    726          CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX  
    727             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout ) 
     729!         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX  
     730!            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout ) 
    728731#if defined key_diahth 
    729732         CALL histdef( nid_T, "sothedep", "Thermocline Depth"                  , "m"      ,   & ! hth 
     
    891894         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping 
    892895      ENDIF 
    893       zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1) 
    894       CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ??? 
     896!      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1) 
     897!      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ??? 
    895898 
    896899#if defined key_diahth 
     
    907910      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress 
    908911 
    909       CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current 
     912      IF( ln_zad_Aimp ) THEN 
     913         CALL histwrite( nid_W, "vovecrtz", it, wn + wi     , ndim_T, ndex_T )    ! vert. current 
     914      ELSE 
     915         CALL histwrite( nid_W, "vovecrtz", it, wn          , ndim_T, ndex_T )    ! vert. current 
     916      ENDIF 
    910917      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef. 
    911918      CALL histwrite( nid_W, "votkeavm", it, avm            , ndim_T, ndex_T )    ! T vert. eddy visc. coef. 
     
    969976      CALL iom_rstput( 0, 0, inum, 'vozocrtx', un                )    ! now i-velocity 
    970977      CALL iom_rstput( 0, 0, inum, 'vomecrty', vn                )    ! now j-velocity 
    971       CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn                )    ! now k-velocity 
     978      IF( ln_zad_Aimp ) THEN 
     979         CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn + wi        )    ! now k-velocity 
     980      ELSE 
     981         CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn             )    ! now k-velocity 
     982      ENDIF 
    972983      IF( ALLOCATED(ahtu) ) THEN 
    973984         CALL iom_rstput( 0, 0, inum,  'ahtu', ahtu              )    ! aht at u-point 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DOM/dommsk.F90

    r11348 r11586  
    100100         &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, & 
    101101         &             cn_ice, nn_ice_dta,                                     & 
    102          &             rn_ice_tem, rn_ice_sal, rn_ice_age,                     & 
    103102         &             ln_vol, nn_volctl, nn_rimwidth 
    104103      !!--------------------------------------------------------------------- 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DOM/domvvl.F90

    r11348 r11586  
    327327      END DO 
    328328      ! 
    329       IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate ! 
    330          !                                                            ! ------baroclinic part------ ! 
     329      IF( (ln_vvl_ztilde .OR. ln_vvl_layer) .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate ! 
     330         !                                                               ! ------baroclinic part------ ! 
    331331         ! I - initialization 
    332332         ! ================== 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DOM/domwri.F90

    r10425 r11586  
    162162      CALL iom_rstput( 0, 0, inum, 'isfdraft', zprt, ktype = jp_r8 )   !    ! nb of ocean T-points 
    163163      !                                                         ! vertical mesh 
    164       CALL iom_rstput( 0, 0, inum, 'e3t_0', e3t_0, ktype = jp_r8  )    !    ! scale factors 
    165       CALL iom_rstput( 0, 0, inum, 'e3u_0', e3u_0, ktype = jp_r8  ) 
    166       CALL iom_rstput( 0, 0, inum, 'e3v_0', e3v_0, ktype = jp_r8  ) 
    167       CALL iom_rstput( 0, 0, inum, 'e3w_0', e3w_0, ktype = jp_r8  ) 
     164      CALL iom_rstput( 0, 0, inum, 'e3t_1d', e3t_1d, ktype = jp_r8  )    !    ! scale factors 
     165      CALL iom_rstput( 0, 0, inum, 'e3w_1d', e3w_1d, ktype = jp_r8  ) 
     166       
     167      CALL iom_rstput( 0, 0, inum, 'e3t_0' , e3t_0 , ktype = jp_r8  ) 
     168      CALL iom_rstput( 0, 0, inum, 'e3u_0' , e3u_0 , ktype = jp_r8  ) 
     169      CALL iom_rstput( 0, 0, inum, 'e3v_0' , e3v_0 , ktype = jp_r8  ) 
     170      CALL iom_rstput( 0, 0, inum, 'e3f_0' , e3f_0 , ktype = jp_r8  ) 
     171      CALL iom_rstput( 0, 0, inum, 'e3w_0' , e3w_0 , ktype = jp_r8  ) 
     172      CALL iom_rstput( 0, 0, inum, 'e3uw_0', e3uw_0, ktype = jp_r8  ) 
     173      CALL iom_rstput( 0, 0, inum, 'e3vw_0', e3vw_0, ktype = jp_r8  ) 
    168174      ! 
    169175      CALL iom_rstput( 0, 0, inum, 'gdept_1d' , gdept_1d , ktype = jp_r8 )  ! stretched system 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DYN/dynhpg.F90

    r11348 r11586  
    3737   USE trd_oce         ! trends: ocean variables 
    3838   USE trddyn          ! trend manager: dynamics 
    39 !jc   USE zpshde          ! partial step: hor. derivative     (zps_hde routine) 
     39   USE zpshde          ! partial step: hor. derivative     (zps_hde routine) 
    4040   ! 
    4141   USE in_out_manager  ! I/O manager 
     
    338338      REAL(wp) ::   zcoef0, zcoef1, zcoef2, zcoef3   ! temporary scalars 
    339339      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zhpi, zhpj 
     340      REAL(wp), DIMENSION(jpi,jpj) :: zgtsu, zgtsv, zgru, zgrv 
    340341      !!---------------------------------------------------------------------- 
    341342      ! 
     
    346347      ENDIF 
    347348 
    348       ! Partial steps: bottom before horizontal gradient of t, s, rd at the last ocean level 
    349 !jc      CALL zps_hde    ( kt, jpts, tsn, gtsu, gtsv, rhd, gru , grv ) 
     349      ! Partial steps: Compute NOW horizontal gradient of t, s, rd at the last ocean level 
     350      CALL zps_hde( kt, jpts, tsn, zgtsu, zgtsv, rhd, zgru , zgrv ) 
    350351 
    351352      ! Local constant initialization 
     
    385386      END DO 
    386387 
    387       ! partial steps correction at the last level  (use gru & grv computed in zpshde.F90) 
     388      ! partial steps correction at the last level  (use zgru & zgrv computed in zpshde.F90) 
    388389      DO jj = 2, jpjm1 
    389390         DO ji = 2, jpim1 
     
    395396               ua  (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku)         ! subtract old value 
    396397               zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1)                   &   ! compute the new one 
    397                   &            + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) * r1_e1u(ji,jj) 
     398                  &            + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + zgru(ji,jj) ) * r1_e1u(ji,jj) 
    398399               ua  (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku)         ! add the new one to the general momentum trend 
    399400            ENDIF 
     
    401402               va  (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv)         ! subtract old value 
    402403               zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1)                   &   ! compute the new one 
    403                   &            + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) * r1_e2v(ji,jj) 
     404                  &            + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + zgrv(ji,jj) ) * r1_e2v(ji,jj) 
    404405               va  (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv)         ! add the new one to the general momentum trend 
    405406            ENDIF 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DYN/sshwzv.F90

    r11413 r11586  
    99   !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module 
    1010   !!            3.3  !  2011-10  (M. Leclair) split former ssh_wzv routine and remove all vvl related work 
     11   !!            4.0  !  2018-12  (A. Coward) add mixed implicit/explicit advection 
    1112   !!---------------------------------------------------------------------- 
    1213 
     
    279280      !!            :   wi      : now vertical velocity (for implicit treatment) 
    280281      !! 
    281       !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
     282      !! Reference  : Shchepetkin, A. F. (2015): An adaptive, Courant-number-dependent 
     283      !!              implicit scheme for vertical advection in oceanic modeling.  
     284      !!              Ocean Modelling, 91, 38-69. 
    282285      !!---------------------------------------------------------------------- 
    283286      INTEGER, INTENT(in) ::   kt   ! time step 
     
    286289      REAL(wp)             ::   zCu, zcff, z1_e3t                     ! local scalars 
    287290      REAL(wp) , PARAMETER ::   Cu_min = 0.15_wp                      ! local parameters 
    288       REAL(wp) , PARAMETER ::   Cu_max = 0.27                         ! local parameters 
     291      REAL(wp) , PARAMETER ::   Cu_max = 0.30_wp                      ! local parameters 
    289292      REAL(wp) , PARAMETER ::   Cu_cut = 2._wp*Cu_max - Cu_min        ! local parameters 
    290293      REAL(wp) , PARAMETER ::   Fcu    = 4._wp*Cu_max*(Cu_max-Cu_min) ! local parameters 
     
    300303      ENDIF 
    301304      ! 
    302       ! 
    303       DO jk = 1, jpkm1            ! calculate Courant numbers 
    304          DO jj = 2, jpjm1 
    305             DO ji = 2, fs_jpim1   ! vector opt. 
    306                z1_e3t = 1._wp / e3t_n(ji,jj,jk) 
    307                Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) )   &  ! 2*rdt and not r2dt (for restartability) 
    308                   &                             + ( MAX( e2u(ji  ,jj)*e3u_n(ji  ,jj,jk)*un(ji  ,jj,jk), 0._wp ) -   & 
    309                   &                                 MIN( e2u(ji-1,jj)*e3u_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) )   & 
    310                   &                               * r1_e1e2t(ji,jj)                                                 & 
    311                   &                             + ( MAX( e1v(ji,jj  )*e3v_n(ji,jj  ,jk)*vn(ji,jj  ,jk), 0._wp ) -   & 
    312                   &                                 MIN( e1v(ji,jj-1)*e3v_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) )   & 
    313                   &                               * r1_e1e2t(ji,jj)                                                 & 
    314                   &                             ) * z1_e3t 
     305      ! Calculate Courant numbers 
     306      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 
     307         DO jk = 1, jpkm1 
     308            DO jj = 2, jpjm1 
     309               DO ji = 2, fs_jpim1   ! vector opt. 
     310                  z1_e3t = 1._wp / e3t_n(ji,jj,jk) 
     311                  ! 2*rdt and not r2dt (for restartability) 
     312                  Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) )                       &   
     313                     &                             + ( MAX( e2u(ji  ,jj)*e3u_n(ji  ,jj,jk)*un(ji  ,jj,jk) + un_td(ji  ,jj,jk), 0._wp ) -   & 
     314                     &                                 MIN( e2u(ji-1,jj)*e3u_n(ji-1,jj,jk)*un(ji-1,jj,jk) + un_td(ji-1,jj,jk), 0._wp ) )   & 
     315                     &                               * r1_e1e2t(ji,jj)                                                                     & 
     316                     &                             + ( MAX( e1v(ji,jj  )*e3v_n(ji,jj  ,jk)*vn(ji,jj  ,jk) + vn_td(ji,jj  ,jk), 0._wp ) -   & 
     317                     &                                 MIN( e1v(ji,jj-1)*e3v_n(ji,jj-1,jk)*vn(ji,jj-1,jk) + vn_td(ji,jj-1,jk), 0._wp ) )   & 
     318                     &                               * r1_e1e2t(ji,jj)                                                                     & 
     319                     &                             ) * z1_e3t 
     320               END DO 
    315321            END DO 
    316322         END DO 
    317       END DO 
     323      ELSE 
     324         DO jk = 1, jpkm1 
     325            DO jj = 2, jpjm1 
     326               DO ji = 2, fs_jpim1   ! vector opt. 
     327                  z1_e3t = 1._wp / e3t_n(ji,jj,jk) 
     328                  ! 2*rdt and not r2dt (for restartability) 
     329                  Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) )   &  
     330                     &                             + ( MAX( e2u(ji  ,jj)*e3u_n(ji  ,jj,jk)*un(ji  ,jj,jk), 0._wp ) -   & 
     331                     &                                 MIN( e2u(ji-1,jj)*e3u_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) )   & 
     332                     &                               * r1_e1e2t(ji,jj)                                                 & 
     333                     &                             + ( MAX( e1v(ji,jj  )*e3v_n(ji,jj  ,jk)*vn(ji,jj  ,jk), 0._wp ) -   & 
     334                     &                                 MIN( e1v(ji,jj-1)*e3v_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) )   & 
     335                     &                               * r1_e1e2t(ji,jj)                                                 & 
     336                     &                             ) * z1_e3t 
     337               END DO 
     338            END DO 
     339         END DO 
     340      ENDIF 
    318341      CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1. ) 
    319342      ! 
     
    346369                  wn(ji,jj,jk) = ( 1._wp - zcff ) * wn(ji,jj,jk) 
    347370                  ! 
    348                   Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient 
     371                  Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient below and in stp_ctl 
    349372               END DO 
    350373            END DO 
     
    353376      ELSE 
    354377         ! Fully explicit everywhere 
    355          Cu_adv(:,:,:) = 0._wp                          ! Reuse array to output coefficient 
     378         Cu_adv(:,:,:) = 0._wp                          ! Reuse array to output coefficient below and in stp_ctl 
    356379         wi    (:,:,:) = 0._wp 
    357380      ENDIF 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/ICB/icbrst.F90

    r10425 r11586  
    131131      CALL iom_get( ncid, jpdom_unknown, 'kount' , zdata(:) ) 
    132132      num_bergs(:) = INT(zdata(:)) 
    133       ! Close file 
    134       CALL iom_close( ncid ) 
    135133      ! 
    136134 
     
    146144      IF( lwp )   WRITE(numout,'(a,i5,a,i5,a)') 'icebergs, icb_rst_read: there were',ibergs_in_file,   & 
    147145         &                                    ' bergs in the restart file and', jn,' bergs have been read' 
     146      ! Close file 
     147      CALL iom_close( ncid ) 
    148148      ! 
    149149      ! Confirm that all areas have a suitable base for assigning new iceberg 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/ICB/icbstp.F90

    r10570 r11586  
    8686      !                                   !* write out time 
    8787      ll_verbose = .FALSE. 
    88       IF( nn_verbose_write > 0 .AND. MOD( kt-1 , nn_verbose_write ) == 0 )   ll_verbose = ( nn_verbose_level >= 0 ) 
     88      IF( nn_verbose_write > 0 .AND. MOD( kt-1 , nn_verbose_write ) == 0 )   ll_verbose = ( nn_verbose_level > 0 ) 
    8989      ! 
    9090      IF( ll_verbose )   WRITE(numicb,9100) nktberg, ndastp, nsec_day 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/IOM/iom.F90

    r11413 r11586  
    3030#if defined key_iomput 
    3131   USE sbc_oce  , ONLY :   nn_fsbc, ght_abl, ghw_abl, e3t_abl, e3w_abl, jpka, jpkam1 
    32    USE trc_oce  , ONLY :   nn_dttrc                  ! frequency of step on passive tracers 
    33    USE icb_oce  , ONLY :   nclasses, class_num       ! iceberg classes 
     32   USE trc_oce  , ONLY :   nn_dttrc        !  !: frequency of step on passive tracers 
     33   USE icb_oce  , ONLY :   nclasses, class_num       !  !: iceberg classes 
    3434#if defined key_si3 
    3535   USE ice      , ONLY :   jpl 
     
    230230          za_bnds(2,:) = ght_abl(2:jpka  ) + e3w_abl(2:jpka) 
    231231          CALL iom_set_axis_attr( "ghw_abl", bounds=za_bnds ) 
     232 
    232233          CALL iom_set_axis_attr( "nfloat", (/ (REAL(ji,wp), ji=1,jpnfl) /) ) 
    233234# if defined key_si3 
     
    19601961      INTEGER :: icnr, jcnr                             ! Offset such that the vertex coordinate (i+icnr,j+jcnr) 
    19611962      !                                                 ! represents the bottom-left corner of cell (i,j) 
    1962       REAL(wp), DIMENSION(4,jpi,jpj,2) ::   z_bnds      ! Lat/lon coordinates of the vertices of cell (i,j) 
    1963       REAL(wp), DIMENSION(  jpi,jpj  ) ::   z_fld       ! Working array to determine where to rotate cells 
    1964       REAL(wp), DIMENSION(4        ,2) ::   z_rot       ! Lat/lon working array for rotation of cells 
    1965       !!---------------------------------------------------------------------- 
     1963      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: z_bnds      ! Lat/lon coordinates of the vertices of cell (i,j) 
     1964      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: z_fld       ! Working array to determine where to rotate cells 
     1965      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: z_rot       ! Lat/lon working array for rotation of cells 
     1966      !!---------------------------------------------------------------------- 
     1967      ! 
     1968      ALLOCATE( z_bnds(4,jpi,jpj,2), z_fld(jpi,jpj), z_rot(4,2)  ) 
    19661969      ! 
    19671970      ! Offset of coordinate representing bottom-left corner 
     
    20432046      CALL iom_set_domain_attr("grid_"//cdgrd, bounds_lat = RESHAPE(z_bnds(:,nldi:nlei,nldj:nlej,1),(/ 4,ni*nj /)),           & 
    20442047          &                                    bounds_lon = RESHAPE(z_bnds(:,nldi:nlei,nldj:nlej,2),(/ 4,ni*nj /)), nvertex=4 ) 
     2048      ! 
     2049      DEALLOCATE( z_bnds, z_fld, z_rot )  
    20452050      ! 
    20462051   END SUBROUTINE set_grid_bounds 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/IOM/iom_nf90.F90

    r11223 r11586  
    196196      CHARACTER(len=*)     , INTENT(in   )           ::   cdvar    ! name of the variable 
    197197      INTEGER              , INTENT(in   )           ::   kiv   !  
    198       INTEGER, DIMENSION(:), INTENT(  out), OPTIONAL ::   kdimsz   ! size of the dimensions 
    199       INTEGER,               INTENT(  out), OPTIONAL ::   kndims   ! size of the dimensions 
     198      INTEGER, DIMENSION(:), INTENT(  out), OPTIONAL ::   kdimsz   ! size of each dimension 
     199      INTEGER              , INTENT(  out), OPTIONAL ::   kndims   ! number of dimensions 
    200200      LOGICAL              , INTENT(  out), OPTIONAL ::   lduld    ! true if the last dimension is unlimited (time) 
    201201      ! 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/LBC/mppini.F90

    r11413 r11586  
    167167           &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, & 
    168168           &             cn_ice, nn_ice_dta,                                     & 
    169            &             rn_ice_tem, rn_ice_sal, rn_ice_age,                     & 
    170169           &             ln_vol, nn_volctl, nn_rimwidth 
    171170      NAMELIST/nammpp/ jpni, jpnj, ln_nnogather, ln_listonly 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcblk.F90

    r11360 r11586  
    121121   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   t_zu, q_zu                 ! air temp. and spec. hum. at wind speed height (L15 bulk scheme) 
    122122 
    123    INTEGER , PUBLIC, PARAMETER ::   jpfld   =11   !: maximum number of files to read 
     123   !INTEGER , PUBLIC, PARAMETER ::   jpfld   =11   !: maximum number of files to read 
     124   INTEGER , PUBLIC            ::   jpfld         !: maximum number of files to read 
    124125   INTEGER , PUBLIC, PARAMETER ::   jp_wndi = 1   !: index of 10m wind velocity (i-component) (m/s)    at T-point 
    125126   INTEGER , PUBLIC, PARAMETER ::   jp_wndj = 2   !: index of 10m wind velocity (j-component) (m/s)    at T-point 
     
    171172      !! 
    172173      CHARACTER(len=100)            ::   cn_dir                ! Root directory for location of atmospheric forcing files 
    173       TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                 ! array of namelist informations on the fields to read 
     174      !TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                 ! array of namelist informations on the fields to read 
     175      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i                 ! array of namelist informations on the fields to read 
    174176      TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_humi, sn_qsr       ! informations about the fields to be read 
    175177      TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow      !       "                        " 
     
    216218      !                                   !* set the bulk structure 
    217219      !                                      !- store namelist information in an array 
     220      IF( ln_blk ) jpfld = 9 
     221      IF( ln_abl ) jpfld = 11 
     222      ALLOCATE( slf_i(jpfld) ) 
     223      ! 
    218224      slf_i(jp_wndi) = sn_wndi   ;   slf_i(jp_wndj) = sn_wndj 
    219225      slf_i(jp_qsr ) = sn_qsr    ;   slf_i(jp_qlw ) = sn_qlw 
     
    221227      slf_i(jp_prec) = sn_prec   ;   slf_i(jp_snow) = sn_snow 
    222228      slf_i(jp_slp ) = sn_slp 
    223       slf_i(jp_hpgi) = sn_hpgi   ;   slf_i(jp_hpgj) = sn_hpgj   
     229      IF( ln_abl ) THEN 
     230        slf_i(jp_hpgi) = sn_hpgi   ;   slf_i(jp_hpgj) = sn_hpgj   
     231      END IF 
    224232      ! 
    225233      !                                      !- allocate the bulk structure 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/TRA/traadv_fct.F90

    r10425 r11586  
    2121   USE diaar5         ! AR5 diagnostics 
    2222   USE phycst  , ONLY : rau0_rcp 
     23   USE zdf_oce , ONLY : ln_zad_Aimp 
    2324   ! 
    2425   USE in_out_manager ! I/O manager 
     
    8687      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw 
    8788      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdx, ztrdy, ztrdz, zptry 
     89      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   zwinf, zwdia, zwsup 
     90      LOGICAL  ::   ll_zAimp                                 ! flag to apply adaptive implicit vertical advection 
    8891      !!---------------------------------------------------------------------- 
    8992      ! 
     
    97100      l_hst = .FALSE. 
    98101      l_ptr = .FALSE. 
     102      ll_zAimp = .FALSE. 
    99103      IF( ( cdtype =='TRA' .AND. l_trdtra  ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) )       l_trd = .TRUE. 
    100104      IF(   cdtype =='TRA' .AND. ln_diaptr )                                                l_ptr = .TRUE.  
     
    116120      ! 
    117121      zwi(:,:,:) = 0._wp         
     122      ! 
     123      ! If adaptive vertical advection, check if it is needed on this PE at this time 
     124      IF( ln_zad_Aimp ) THEN 
     125         IF( MAXVAL( ABS( wi(:,:,:) ) ) > 0._wp ) ll_zAimp = .TRUE. 
     126      END IF 
     127      ! If active adaptive vertical advection, build tridiagonal matrix 
     128      IF( ll_zAimp ) THEN 
     129         ALLOCATE(zwdia(jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk)) 
     130         DO jk = 1, jpkm1 
     131            DO jj = 2, jpjm1 
     132               DO ji = fs_2, fs_jpim1   ! vector opt. (ensure same order of calculation as below if wi=0.) 
     133                  zwdia(ji,jj,jk) =  1._wp + p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) / e3t_a(ji,jj,jk) 
     134                  zwinf(ji,jj,jk) =  p2dt * MIN( wi(ji,jj,jk  ) , 0._wp ) / e3t_a(ji,jj,jk) 
     135                  zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t_a(ji,jj,jk) 
     136               END DO 
     137            END DO 
     138         END DO 
     139      END IF 
    118140      ! 
    119141      DO jn = 1, kjpt            !==  loop over the tracers  ==! 
     
    169191            END DO 
    170192         END DO 
     193          
     194         IF ( ll_zAimp ) THEN 
     195            CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 ) 
     196            ! 
     197            ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ; 
     198            DO jk = 2, jpkm1        ! Interior value ( multiplied by wmask) 
     199               DO jj = 2, jpjm1 
     200                  DO ji = fs_2, fs_jpim1   ! vector opt.   
     201                     zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 
     202                     zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 
     203                     ztw(ji,jj,jk) =  0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk) 
     204                     zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! update vertical fluxes 
     205                  END DO 
     206               END DO 
     207            END DO 
     208            DO jk = 1, jpkm1 
     209               DO jj = 2, jpjm1 
     210                  DO ji = fs_2, fs_jpim1   ! vector opt.   
     211                     pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) & 
     212                        &                                  * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     213                  END DO 
     214               END DO 
     215            END DO 
     216            ! 
     217         END IF 
    171218         !                 
    172219         IF( l_trd .OR. l_hst )  THEN             ! trend diagnostics (contribution of upstream fluxes) 
     
    277324            zwz(:,:,1) = 0._wp   ! only ocean surface as interior zwz values have been w-masked 
    278325         ENDIF 
     326         !          
     327         IF ( ll_zAimp ) THEN 
     328            DO jk = 1, jpkm1     !* trend and after field with monotonic scheme 
     329               DO jj = 2, jpjm1 
     330                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     331                     !                             ! total intermediate advective trends 
     332                     ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
     333                        &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
     334                        &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj) 
     335                     ztw(ji,jj,jk)  = zwi(ji,jj,jk) + p2dt * ztra / e3t_a(ji,jj,jk) * tmask(ji,jj,jk) 
     336                  END DO 
     337               END DO 
     338            END DO 
     339            ! 
     340            CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 ) 
     341            ! 
     342            DO jk = 2, jpkm1        ! Interior value ( multiplied by wmask) 
     343               DO jj = 2, jpjm1 
     344                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     345                     zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 
     346                     zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 
     347                     zwz(ji,jj,jk) =  zwz(ji,jj,jk) + 0.5 * e1e2t(ji,jj) * ( zfp_wk * ztw(ji,jj,jk) + zfm_wk * ztw(ji,jj,jk-1) ) * wmask(ji,jj,jk) 
     348                  END DO 
     349               END DO 
     350            END DO 
     351         END IF 
    279352         ! 
    280353         CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1., zwx, 'U', -1. , zwy, 'V', -1.,  zwz, 'W',  1. ) 
     
    289362            DO jj = 2, jpjm1 
    290363               DO ji = fs_2, fs_jpim1   ! vector opt.   
    291                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
    292                      &                                   + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
    293                      &                                   + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) & 
    294                      &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    295                END DO 
    296             END DO 
    297          END DO 
     364                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
     365                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
     366                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj) 
     367                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra / e3t_n(ji,jj,jk) 
     368                  zwi(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t_a(ji,jj,jk) * tmask(ji,jj,jk) 
     369               END DO 
     370            END DO 
     371         END DO 
     372         ! 
     373         IF ( ll_zAimp ) THEN 
     374            ! 
     375            ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp 
     376            DO jk = 2, jpkm1        ! Interior value ( multiplied by wmask) 
     377               DO jj = 2, jpjm1 
     378                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     379                     zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 
     380                     zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 
     381                     ztw(ji,jj,jk) = - 0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk) 
     382                     zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! Update vertical fluxes for trend diagnostic 
     383                  END DO 
     384               END DO 
     385            END DO 
     386            DO jk = 1, jpkm1 
     387               DO jj = 2, jpjm1 
     388                  DO ji = fs_2, fs_jpim1   ! vector opt.   
     389                     pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) & 
     390                        &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     391                  END DO 
     392               END DO 
     393            END DO 
     394         END IF          
    298395         ! 
    299396         IF( l_trd .OR. l_hst ) THEN   ! trend diagnostics // heat/salt transport 
     
    318415      END DO                     ! end of tracer loop 
    319416      ! 
     417      IF ( ll_zAimp ) THEN 
     418         DEALLOCATE( zwdia, zwinf, zwsup ) 
     419      ENDIF 
    320420      IF( l_trd .OR. l_hst ) THEN  
    321421         DEALLOCATE( ztrdx, ztrdy, ztrdz ) 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/TRA/traldf_iso.F90

    r10068 r11586  
    289289         !!---------------------------------------------------------------------- 
    290290         ! 
    291          ztfw(1,:,:) = 0._wp     ;     ztfw(jpi,:,:) = 0._wp 
     291         ztfw(fs_2:1,:,:) = 0._wp     ;     ztfw(jpi:fs_jpim1,:,:) = 0._wp   ! avoid to potentially manipulate NaN values 
    292292         ! 
    293293         ! Vertical fluxes 
     
    323323         IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz 
    324324            DO jk = 2, jpkm1        
    325                DO jj = 1, jpjm1 
     325               DO jj = 2, jpjm1 
    326326                  DO ji = fs_2, fs_jpim1 
    327327                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)   & 
     
    336336            CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2 
    337337               DO jk = 2, jpkm1  
    338                   DO jj = 1, jpjm1 
     338                  DO jj = 2, jpjm1 
    339339                     DO ji = fs_2, fs_jpim1 
    340340                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk)    & 
     
    346346            CASE(  2  )                         ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb  and ptbb gradients, resp. 
    347347               DO jk = 2, jpkm1  
    348                   DO jj = 1, jpjm1 
     348                  DO jj = 2, jpjm1 
    349349                     DO ji = fs_2, fs_jpim1 
    350350                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)                      & 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/step.F90

    r11413 r11586  
    165165                            CALL eos    ( tsn, rhd, rhop, gdept_n(:,:,:) )  ! now in situ density for hpg computation 
    166166                             
    167 !!jc: fs simplification 
    168 !!jc: lines below are useless if ln_linssh=F. Keep them here (which maintains a bug if ln_linssh=T and ln_zps=T, cf ticket #1636)  
    169 !!                                         but ensures reproductible results 
    170 !!                                         with previous versions using split-explicit free surface           
    171             IF( ln_zps .AND. .NOT. ln_isfcav )                               & 
    172                &            CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,   &  ! Partial steps: before horizontal gradient 
    173                &                                          rhd, gru , grv     )  ! of t, s, rd at the last ocean level 
    174             IF( ln_zps .AND.       ln_isfcav )                                          & 
    175                &            CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF) 
    176                &                                          rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level 
    177 !!jc: fs simplification 
    178167                             
    179168                         ua(:,:,:) = 0._wp            ! set dynamics trends to zero 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/stpctl.F90

    r10570 r11586  
    9696            IF( ln_zad_Aimp ) THEN 
    9797               istatus = NF90_DEF_VAR( idrun,   'abs_wi_max', NF90_DOUBLE, (/ idtime /), idw1  ) 
    98                istatus = NF90_DEF_VAR( idrun,       'Cu_max', NF90_DOUBLE, (/ idtime /), idc1  ) 
     98               istatus = NF90_DEF_VAR( idrun,       'Cf_max', NF90_DOUBLE, (/ idtime /), idc1  ) 
    9999            ENDIF 
    100100            istatus = NF90_ENDDEF(idrun) 
     
    123123      IF( ln_zad_Aimp ) THEN 
    124124         zmax(8) = MAXVAL(  ABS( wi(:,:,:) ) , mask = wmask(:,:,:) == 1._wp ) ! implicit vertical vel. max 
    125          zmax(9) = MAXVAL(   Cu_adv(:,:,:)   , mask = tmask(:,:,:) == 1._wp ) !       cell Courant no. max 
     125         zmax(9) = MAXVAL(   Cu_adv(:,:,:)   , mask = tmask(:,:,:) == 1._wp ) ! partitioning coeff. max 
    126126      ENDIF 
    127127      ! 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/tests/BENCH/EXPREF/best_jpni_jpnj_eorca025

    r9856 r11586  
    624624    nb_cores 36180 ( 180 x 201 ), nb_points 80 ( 10 x 8 ) 
    625625    nb_cores 38560 ( 160 x 241 ), nb_points 77 ( 11 x 7 ) 
     626    nb_cores 41406 ( 206 x 201 ), nb_points 72 ( 9 x 8 ) 
     627    nb_cores 43380 ( 180 x 241 ), nb_points 70 ( 10 x 7 ) 
     628    nb_cores 48240 ( 240 x 201 ), nb_points 64 ( 8 x 8 ) 
     629    nb_cores 49646 ( 206 x 241 ), nb_points 63 ( 9 x 7 ) 
     630    nb_cores 54360 ( 180 x 302 ), nb_points 60 ( 10 x 6 ) 
     631    nb_cores 54360 ( 360 x 151 ), nb_points 60 ( 6 x 10 ) 
     632    nb_cores 57840 ( 240 x 241 ), nb_points 56 ( 8 x 7 ) 
     633    nb_cores 62212 ( 206 x 302 ), nb_points 54 ( 9 x 6 ) 
     634    nb_cores 69408 ( 288 x 241 ), nb_points 49 ( 7 x 7 ) 
     635    nb_cores 72360 ( 360 x 201 ), nb_points 48 ( 6 x 8 ) 
     636    nb_cores 82812 ( 206 x 402 ), nb_points 45 ( 9 x 5 ) 
     637    nb_cores 86760 ( 360 x 241 ), nb_points 42 ( 6 x 7 ) 
     638    nb_cores 96480 ( 240 x 402 ), nb_points 40 ( 8 x 5 ) 
     639    nb_cores 96480 ( 480 x 201 ), nb_points 40 ( 5 x 8 ) 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/tests/BENCH/EXPREF/best_jpni_jpnj_eorca12

    r9856 r11586  
    13001300    nb_cores 49852 ( 206 x 242 ), nb_points 345 ( 23 x 15 ) 
    13011301    nb_cores 49950 ( 270 x 185 ), nb_points 342 ( 18 x 19 ) 
     1302    nb_cores 50400 ( 288 x 175 ), nb_points 340 ( 17 x 20 ) 
     1303    nb_cores 50400 ( 240 x 210 ), nb_points 340 ( 20 x 17 ) 
     1304    nb_cores 51294 ( 309 x 166 ), nb_points 336 ( 16 x 21 ) 
     1305    nb_cores 52272 ( 216 x 242 ), nb_points 330 ( 22 x 15 ) 
     1306    nb_cores 53190 ( 270 x 197 ), nb_points 324 ( 18 x 18 ) 
     1307    nb_cores 53280 ( 288 x 185 ), nb_points 323 ( 17 x 19 ) 
     1308    nb_cores 54000 ( 240 x 225 ), nb_points 320 ( 20 x 16 ) 
     1309    nb_cores 55176 ( 228 x 242 ), nb_points 315 ( 21 x 15 ) 
     1310    nb_cores 56199 ( 393 x 143 ), nb_points 312 ( 13 x 24 ) 
     1311    nb_cores 56700 ( 270 x 210 ), nb_points 306 ( 18 x 17 ) 
     1312    nb_cores 57165 ( 309 x 185 ), nb_points 304 ( 16 x 19 ) 
     1313    nb_cores 58080 ( 240 x 242 ), nb_points 300 ( 20 x 15 ) 
     1314    nb_cores 58916 ( 206 x 286 ), nb_points 299 ( 23 x 13 ) 
     1315    nb_cores 59760 ( 360 x 166 ), nb_points 294 ( 14 x 21 ) 
     1316    nb_cores 60480 ( 288 x 210 ), nb_points 289 ( 17 x 17 ) 
     1317    nb_cores 60750 ( 270 x 225 ), nb_points 288 ( 18 x 16 ) 
     1318    nb_cores 61605 ( 333 x 185 ), nb_points 285 ( 15 x 19 ) 
     1319    nb_cores 63000 ( 360 x 175 ), nb_points 280 ( 14 x 20 ) 
     1320    nb_cores 64800 ( 288 x 225 ), nb_points 272 ( 17 x 16 ) 
     1321    nb_cores 65340 ( 270 x 242 ), nb_points 270 ( 18 x 15 ) 
     1322    nb_cores 66600 ( 360 x 185 ), nb_points 266 ( 14 x 19 ) 
     1323    nb_cores 68040 ( 216 x 315 ), nb_points 264 ( 22 x 12 ) 
     1324    nb_cores 68640 ( 240 x 286 ), nb_points 260 ( 20 x 13 ) 
     1325    nb_cores 69525 ( 309 x 225 ), nb_points 256 ( 16 x 16 ) 
     1326    nb_cores 69696 ( 288 x 242 ), nb_points 255 ( 17 x 15 ) 
     1327    nb_cores 70920 ( 360 x 197 ), nb_points 252 ( 14 x 18 ) 
     1328    nb_cores 72705 ( 393 x 185 ), nb_points 247 ( 13 x 19 ) 
     1329    nb_cores 74778 ( 309 x 242 ), nb_points 240 ( 16 x 15 ) 
     1330    nb_cores 75600 ( 360 x 210 ), nb_points 238 ( 14 x 17 ) 
     1331    nb_cores 77220 ( 270 x 286 ), nb_points 234 ( 18 x 13 ) 
     1332    nb_cores 79680 ( 480 x 166 ), nb_points 231 ( 11 x 21 ) 
     1333    nb_cores 79920 ( 432 x 185 ), nb_points 228 ( 12 x 19 ) 
     1334    nb_cores 80586 ( 333 x 242 ), nb_points 225 ( 15 x 15 ) 
     1335    nb_cores 81000 ( 360 x 225 ), nb_points 224 ( 14 x 16 ) 
     1336    nb_cores 82368 ( 288 x 286 ), nb_points 221 ( 17 x 13 ) 
     1337    nb_cores 84000 ( 480 x 175 ), nb_points 220 ( 11 x 20 ) 
     1338    nb_cores 84000 ( 240 x 350 ), nb_points 220 ( 20 x 11 ) 
     1339    nb_cores 85050 ( 270 x 315 ), nb_points 216 ( 18 x 12 ) 
     1340    nb_cores 87120 ( 360 x 242 ), nb_points 210 ( 14 x 15 ) 
     1341    nb_cores 88374 ( 309 x 286 ), nb_points 208 ( 16 x 13 ) 
     1342    nb_cores 90720 ( 288 x 315 ), nb_points 204 ( 17 x 12 ) 
     1343    nb_cores 90720 ( 432 x 210 ), nb_points 204 ( 12 x 17 ) 
     1344    nb_cores 94500 ( 270 x 350 ), nb_points 198 ( 18 x 11 ) 
     1345    nb_cores 94680 ( 360 x 263 ), nb_points 196 ( 14 x 14 ) 
     1346    nb_cores 95106 ( 393 x 242 ), nb_points 195 ( 13 x 15 ) 
     1347    nb_cores 97200 ( 432 x 225 ), nb_points 192 ( 12 x 16 ) 
     1348    nb_cores 99900 ( 540 x 185 ), nb_points 190 ( 10 x 19 ) 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/tests/BENCH/EXPREF/namelist_cfg_orca025_like

    r10539 r11586  
    99   nn_it000    =       1   !  first time step 
    1010   nn_itend    =    1000   !  last time step  
    11    nn_stock    =       0   !  frequency of creation of a restart file (modulo referenced to 1) 
    12    nn_write    =       0   !  frequency of write in the output file   (modulo referenced to nn_it000) 
     11   nn_stock    =      -1   !  frequency of creation of a restart file (modulo referenced to 1) 
     12   nn_write    =      -1   !  frequency of write in the output file   (modulo referenced to nn_it000) 
    1313/ 
    1414!----------------------------------------------------------------------- 
     
    3030&namctl        !   Control prints                                       (default: OFF) 
    3131!----------------------------------------------------------------------- 
    32    ln_ctl      = .false.   !  trends control print (expensive!) 
    3332   nn_print    =    0      !  level of print (0 no extra print) 
    3433   ln_timing   = .false.   !  timing by routine write out in timing.output file 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/tests/BENCH/EXPREF/namelist_cfg_orca12_like

    r10539 r11586  
    99   nn_it000    =       1   !  first time step 
    1010   nn_itend    =    1000   !  last time step  
    11    nn_stock    =       0   !  frequency of creation of a restart file (modulo referenced to 1) 
    12    nn_write    =       0   !  frequency of write in the output file   (modulo referenced to nn_it000) 
     11   nn_stock    =      -1   !  frequency of creation of a restart file (modulo referenced to 1) 
     12   nn_write    =      -1   !  frequency of write in the output file   (modulo referenced to nn_it000) 
    1313/ 
    1414!----------------------------------------------------------------------- 
     
    3030&namctl        !   Control prints                                       (default: OFF) 
    3131!----------------------------------------------------------------------- 
    32    ln_ctl      = .false.   !  trends control print (expensive!) 
    3332   nn_print    =    0      !  level of print (0 no extra print) 
    3433   ln_timing   = .false.   !  timing by routine write out in timing.output file 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/tests/BENCH/EXPREF/namelist_cfg_orca1_like

    r10539 r11586  
    99   nn_it000    =       1   !  first time step 
    1010   nn_itend    =    1000   !  last time step  
    11    nn_stock    =       0   !  frequency of creation of a restart file (modulo referenced to 1) 
    12    nn_write    =       0   !  frequency of write in the output file   (modulo referenced to nn_it000) 
     11   nn_stock    =      -1   !  frequency of creation of a restart file (modulo referenced to 1) 
     12   nn_write    =      -1   !  frequency of write in the output file   (modulo referenced to nn_it000) 
    1313/ 
    1414!----------------------------------------------------------------------- 
     
    3030&namctl        !   Control prints                                       (default: OFF) 
    3131!----------------------------------------------------------------------- 
    32    ln_ctl      = .false.   !  trends control print (expensive!) 
    3332   nn_print    =    0      !  level of print (0 no extra print) 
    3433   ln_timing   = .false.   !  timing by routine write out in timing.output file 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/tests/BENCH/EXPREF/namelist_ice_cfg

    r10179 r11586  
    7474&namini         !   Ice initialization 
    7575!------------------------------------------------------------------------------ 
     76   rn_thres_sst     =   0.5           !  max delta temp. above Tfreeze with initial ice = (sst - tfreeze) 
    7677/ 
    7778!------------------------------------------------------------------------------ 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/tests/BENCH/EXPREF/sh_bench

    r9869 r11586  
    99machine=$( hostname | sed -e "s/[0-9]*//g" ) 
    1010case $machine in  
     11    "jean-zay")     ncore_node=40 ;; 
    1112    "beaufixlogin") ncore_node=40 ;; 
    1213    "curie")        ncore_node=16 ;; 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/tests/BENCH/EXPREF/submit_bench

    r9870 r11586  
    66# 
    77set -u 
    8 set -vx 
     8#set -vx 
    99# 
    1010cores=$1 
     
    1515# 
    1616# number of processes for each executable 
    17     nproc_exe1=$( echo $cores | bc ) 
    18   (( nproc = $nproc_exe1 )) 
    19 (( nnode = $nproc / $ncore_node )) 
    20 [ $nnode -lt 1 ] && nnode=1 
     17nproc_exe1=$( echo $cores | bc ) 
     18nproc=$nproc_exe1 
     19nnode=$(( ( $nproc + $ncore_node - 1 )  / $ncore_node )) 
    2120 
    2221nproc5=$( printf "%05d\n" ${nproc_exe1} ) 
     22 
     23case ${resolution} in 
     24    "1") 
     25   if [ $nproc_exe1 -lt 50 ] 
     26   then 
     27       timejob=3600 
     28   elif [ $nproc_exe1 -lt 100 ] 
     29   then 
     30       timejob=1800 
     31   elif [ $nproc_exe1 -lt 200 ] 
     32   then 
     33       timejob=900 
     34   else 
     35       timejob=600 
     36   fi 
     37   ;; 
     38    "025") 
     39   if [ $nproc_exe1 -lt 50 ] 
     40   then 
     41       timejob=15000 
     42   elif [ $nproc_exe1 -lt 100 ] 
     43   then 
     44       timejob=7000 
     45   elif [ $nproc_exe1 -lt 200 ] 
     46   then 
     47       timejob=3600 
     48   elif [ $nproc_exe1 -lt 400 ] 
     49   then 
     50       timejob=2000 
     51   elif [ $nproc_exe1 -lt 800 ] 
     52   then 
     53       timejob=1000 
     54   else 
     55       timejob=600 
     56   fi 
     57   ;; 
     58    "12") 
     59   if [ $nproc_exe1 -lt 200 ] 
     60   then 
     61       timejob=30000 
     62   elif [ $nproc_exe1 -lt 400 ] 
     63   then 
     64       timejob=15000 
     65   elif [ $nproc_exe1 -lt 800 ] 
     66   then 
     67       timejob=20000 
     68   elif [ $nproc_exe1 -lt 1600 ] 
     69   then 
     70       timejob=15000 
     71   elif [ $nproc_exe1 -lt 3200 ] 
     72   then 
     73       timejob=7500 
     74   elif [ $nproc_exe1 -lt 10000 ] 
     75   then 
     76       timejob=5000 
     77   elif [ $nproc_exe1 -lt 20000 ] 
     78   then 
     79       timejob=2500 
     80   else 
     81       timejob=1200 
     82   fi 
     83   ;; 
     84esac 
     85 
    2386 
    2487###################################################################### 
     
    100163if [[ ( "$machine" == "curie" ) || ( "$machine" == "irene" ) ]] 
    101164then 
    102      
    103     case ${resolution} in 
    104    "1") 
    105        if [ $nproc_exe1 -lt 50 ] 
    106        then 
    107       timejob=3600 
    108        elif [ $nproc_exe1 -lt 100 ] 
    109        then 
    110       timejob=1800 
    111        elif [ $nproc_exe1 -lt 200 ] 
    112        then 
    113       timejob=900 
    114        else 
    115       timejob=600 
    116        fi 
    117        ;; 
    118    "025") 
    119        if [ $nproc_exe1 -lt 50 ] 
    120        then 
    121       timejob=15000 
    122        elif [ $nproc_exe1 -lt 100 ] 
    123        then 
    124       timejob=7000 
    125        elif [ $nproc_exe1 -lt 200 ] 
    126        then 
    127       timejob=3600 
    128        elif [ $nproc_exe1 -lt 400 ] 
    129        then 
    130       timejob=2000 
    131        elif [ $nproc_exe1 -lt 800 ] 
    132        then 
    133       timejob=1000 
    134        else 
    135       timejob=600 
    136        fi 
    137        ;; 
    138    "12") 
    139        if [ $nproc_exe1 -lt 200 ] 
    140        then 
    141       timejob=30000 
    142        elif [ $nproc_exe1 -lt 400 ] 
    143        then 
    144       timejob=15000 
    145        elif [ $nproc_exe1 -lt 800 ] 
    146        then 
    147       timejob=8000 
    148        elif [ $nproc_exe1 -lt 1600 ] 
    149        then 
    150       timejob=5000 
    151        elif [ $nproc_exe1 -lt 3200 ] 
    152        then 
    153       timejob=2500 
    154        else 
    155       timejob=1200 
    156        fi 
    157    ;; 
    158     esac 
    159165 
    160166    [ "$machine" == "curie" ] && queuename=standard || queuename=skylake 
    161167     
    162     jobname=$HOME/binrun/jobbench 
     168    EXPjob=../EXP_${resolution}_${nproc5}_${dateref} 
     169    mkdir -p ${EXPjob} 
     170    cd ${EXPjob} 
     171    jobname=jobbench 
    163172    cat > $jobname << EOF 
    164173#!/bin/bash 
     
    176185# 
    177186 
    178 cd \${BRIDGE_MSUB_PWD}/.. 
    179  
    180  
    181 EXPjob=EXP_${resolution}_${nproc5}_${dateref} 
    182 rsync -av --exclude="*eo" EXPREF/ \${EXPjob}/ 
    183 rsync -av EXP00/nemo \${EXPjob}/nemo 
    184 cd \${EXPjob} 
     187cd \${BRIDGE_MSUB_PWD} 
     188 
     189for ff in \${BRIDGE_MSUB_PWD}/../EXPREF/namelist_*cfg \${BRIDGE_MSUB_PWD}/../EXPREF/namelist_*ref \${BRIDGE_MSUB_PWD}/../BLD/bin/nemo.exe 
     190do 
     191    cp \$ff . 
     192done 
    185193 
    186194jpni=${cores/\**/} 
    187195jpnj=${cores/?*\*/} 
    188196 
    189 sed -e "s/jpni *=.*/jpni = \${jpni}/" -e "s/jpnj *=.*/jpnj = \${jpnj}/" namelist_cfg_orca${resolution}_like > namelist_cfg 
    190  
    191 time ccc_mprun -n \${BRIDGE_MSUB_NPROC} ./nemo > jobout_${resolution}_${nproc5}_${dateref} 
     197sed -e "s/jpni *=.*/jpni = \${jpni}/" \ 
     198    -e "s/jpnj *=.*/jpnj = \${jpnj}/"\ 
     199    -e "s/ln_timing *= *.false./ln_timing   =  .true./" \ 
     200     \${BRIDGE_MSUB_PWD}/../EXPREF/namelist_cfg_orca${resolution}_like > namelist_cfg 
     201 
     202time ccc_mprun -n \${BRIDGE_MSUB_NPROC} ./nemo.exe > jobout_${resolution}_${nproc5}_${dateref} 2>&1 
    192203 
    193204EOF 
     
    196207 
    197208fi 
     209 
     210###################################################################### 
     211### Jean-Zay 
     212###################################################################### 
     213 
     214if [ "$machine" == "jean-zay" ] 
     215then 
     216    hh=$( printf "%02d\n" $(( ${timejob} / 3600 )) ) 
     217    mm=$( printf "%02d\n" $(( ( ${timejob} % 3600 ) / 60 )) ) 
     218    ss=$( printf "%02d\n" $(( ( ${timejob} % 3600 ) % 60 )) ) 
     219 
     220    EXPjob=../EXP_${resolution}_${nproc5}_${dateref} 
     221    mkdir -p ${EXPjob} 
     222    cd ${EXPjob} 
     223    jobname=jobbench 
     224    cat > $jobname << EOF 
     225#!/bin/bash 
     226#SBATCH --job-name=Seq        # nom du job 
     227#SBATCH --partition=cpu_gct3   # demande d'allocation sur la partition CPU 
     228#SBATCH --nodes=${nnode}             # nombre de noeuds 
     229#SBATCH --ntasks-per-node=${ncore_node}  # nombre de taches MPI par noeud 
     230#SBATCH --ntasks-per-core=1        # 1 processus MPI par coeur physique (pas d'hyperthreading) 
     231#SBATCH --time=${hh}:${mm}:${ss}       # temps d execution maximum demande (HH:MM:SS) 
     232#SBATCH --output=bench_${resolution}_${nproc5}_%j.eo        # nom du fichier de sortie 
     233#SBATCH --error=bench_${resolution}_${nproc5}_%j.eo         # nom du fichier d'erreur (ici en commun avec la sortie) 
     234#========================================== 
     235set -u 
     236#set -xv 
     237# 
     238#cd \${SLURM_SUBMIT_DIR} 
     239cd \${JOBSCRATCH} 
     240pwd 
     241 
     242for ff in \${SLURM_SUBMIT_DIR}/../EXPREF/namelist_*cfg \${SLURM_SUBMIT_DIR}/../EXPREF/namelist_*ref \${SLURM_SUBMIT_DIR}/../BLD/bin/nemo.exe 
     243do 
     244    cp \$ff . 
     245done 
     246 
     247jpni=${cores/\**/} 
     248jpnj=${cores/?*\*/} 
     249 
     250sed -e "s/jpni *=.*/jpni = \${jpni}/" \ 
     251    -e "s/jpnj *=.*/jpnj = \${jpnj}/" \ 
     252    -e "s/ln_timing *= *.false./ln_timing   =  .true./" \ 
     253    \${SLURM_SUBMIT_DIR}/../EXPREF/namelist_cfg_orca${resolution}_like > namelist_cfg 
     254 
     255ls -l 
     256 
     257echo  
     258echo  
     259echo " =========== start the model ===========" 
     260echo  
     261echo  
     262 
     263time srun --mpi=pmi2 --cpu-bind=cores -K1 -n ${nproc} ./nemo.exe > jobout_${resolution}_${nproc5}_${dateref} 2>&1 
     264 
     265ls -l 
     266 
     267if [ "\$( pwd )" != "\${SLURM_SUBMIT_DIR}" ] 
     268then  
     269    rsync -av namelist_cfg time.step ocean.output jobout_${resolution}_${nproc5}_${dateref} communication_report.txt layout.dat timing.output output.namelist* \${SLURM_SUBMIT_DIR} 
     270fi 
     271 
     272EOF 
     273 
     274    sbatch $jobname 
     275 
     276fi 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/tests/BENCH/MY_SRC/usrdef_istate.F90

    r10179 r11586  
    6464      DO jj = 1, jpj 
    6565         DO ji = 1, jpi 
    66             z2d(ji,jj) = 0.1 * ( 0.5 - REAL( nimpp + ji - 1 + ( njmpp + jj - 2 ) * jpiglo, wp ) / REAL( jpiglo * jpjglo, wp ) ) 
     66            z2d(ji,jj) = 0.1 * ( 0.5 - REAL( mig(ji) + mjg(jj) * jpiglo, wp ) / REAL( jpiglo * jpjglo, wp ) ) 
    6767         ENDDO 
    6868      ENDDO 
     
    7373      DO jk = 1, jpk 
    7474         zfact = REAL(jk-1,wp) / REAL(jpk-1,wp)   ! 0 to 1 to add a basic stratification 
    75          ! temperature choosen to lead to 20% ice 
    76          pts(:,:,jk,jp_tem) = 2._wp - 0.1_wp * zfact + z2d(:,:) * 100._wp ! 2 to 1.9 +/- 5 degG 
    77          WHERE ( pts(:,:,jk,jp_tem) < -1.5_wp ) pts(:,:,jk,jp_tem) = -1.5_wp + z2d(:,:) * 0.2_wp   
     75         ! temperature choosen to lead to ~50% ice at the beginning if rn_thres_sst = 0.5 
     76         pts(:,:,jk,jp_tem) = 20._wp*z2d(:,:) - 1._wp - 0.5_wp * zfact    ! -1 to -1.5 +/-1.0 degG 
    7877         ! salinity:   
    7978         pts(:,:,jk,jp_sal) = 30._wp + 1._wp * zfact + z2d(:,:)           ! 30 to 31 +/- 0.05 psu 
     
    8483      ! 
    8584      CALL lbc_lnk('usrdef_istate', pssh, 'T',  1. )            ! apply boundary conditions 
    86       CALL lbc_lnk( 'usrdef_istate', pts, 'T',  1. )            ! apply boundary conditions 
    87       CALL lbc_lnk(  'usrdef_istate', pu, 'U', -1. )            ! apply boundary conditions 
    88       CALL lbc_lnk(  'usrdef_istate', pv, 'V', -1. )            ! apply boundary conditions 
     85      CALL lbc_lnk('usrdef_istate', pts, 'T',  1. )            ! apply boundary conditions 
     86      CALL lbc_lnk('usrdef_istate',  pu, 'U', -1. )            ! apply boundary conditions 
     87      CALL lbc_lnk('usrdef_istate',  pv, 'V', -1. )            ! apply boundary conditions 
    8988       
    9089   END SUBROUTINE usr_def_istate 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/tests/BENCH/MY_SRC/usrdef_nam.F90

    r10428 r11586  
    2929CONTAINS 
    3030 
    31    SUBROUTINE usr_def_nam( ldtxt, ldnam, cd_cfg, kk_cfg, kpi, kpj, kpk, kperio ) 
     31   SUBROUTINE usr_def_nam( cd_cfg, kk_cfg, kpi, kpj, kpk, kperio ) 
    3232      !!---------------------------------------------------------------------- 
    3333      !!                     ***  ROUTINE dom_nam  *** 
     
    4141      !! ** input   : - namusr_def namelist found in namelist_cfg 
    4242      !!---------------------------------------------------------------------- 
    43       CHARACTER(len=*), DIMENSION(:), INTENT(out) ::   ldtxt, ldnam    ! stored print information 
    4443      CHARACTER(len=*)              , INTENT(out) ::   cd_cfg          ! configuration name 
    4544      INTEGER                       , INTENT(out) ::   kk_cfg          ! configuration resolution 
     
    4847      ! 
    4948      ! 
    50       INTEGER ::   ios, ii      ! Local integer 
     49      INTEGER ::   ios         ! Local integer 
    5150      !                              !!* namusr_def namelist *!! 
    5251      INTEGER ::   nn_isize    ! number of point in i-direction of global(local) domain if >0 (<0)   
     
    5554      INTEGER ::   nn_perio    ! periodicity 
    5655      !                              !!* nammpp namelist *!! 
    57       CHARACTER(len=1) ::   cn_mpi_send 
    58       INTEGER          ::   nn_buffer, jpni, jpnj 
     56      INTEGER          ::   jpni, jpnj 
    5957      LOGICAL          ::   ln_nnogather 
    6058      !! 
    6159      NAMELIST/namusr_def/ nn_isize, nn_jsize, nn_ksize, nn_perio 
    62       NAMELIST/nammpp/ cn_mpi_send, nn_buffer, jpni, jpnj, ln_nnogather 
     60      NAMELIST/nammpp/ jpni, jpnj, ln_nnogather 
    6361      !!----------------------------------------------------------------------      
    6462      ! 
    6563      REWIND( numnam_cfg )          ! Namelist namusr_def (exist in namelist_cfg only) 
    6664      READ  ( numnam_cfg, namusr_def, IOSTAT = ios, ERR = 903 ) 
    67 903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namusr_def in configuration namelist', .TRUE. ) 
    68       WRITE( ldnam(:), namusr_def )       
     65903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namusr_def in configuration namelist' ) 
     66      ! 
     67      IF(lwm)   WRITE( numond, namusr_def )       
    6968      ! 
    7069      cd_cfg = 'BENCH'             ! name & resolution (not used) 
    7170      kk_cfg = 0 
    72  
     71      ! 
    7372      IF( nn_isize < 0 .AND. nn_jsize < 0 ) THEN 
    7473      ! 
    7574         REWIND( numnam_ref )              ! Namelist nammpp in reference namelist: mpi variables 
    7675         READ  ( numnam_ref, nammpp, IOSTAT = ios, ERR = 901) 
    77 901      IF( ios /= 0 )   CALL ctl_nam ( ios , 'nammpp in reference namelist', lwp ) 
     76901      IF( ios /= 0 )   CALL ctl_nam ( ios , 'nammpp in reference namelist' ) 
    7877         ! 
    7978         REWIND( numnam_cfg )              ! Namelist nammpp in configuration namelist: mpi variables 
    8079         READ  ( numnam_cfg, nammpp, IOSTAT = ios, ERR = 902 ) 
    81 902      IF( ios >  0 )   CALL ctl_nam ( ios , 'nammpp in configuration namelist', lwp ) 
     80902      IF( ios >  0 )   CALL ctl_nam ( ios , 'nammpp in configuration namelist' ) 
    8281 
    8382         kpi = ( -nn_isize - 2*nn_hls ) * jpni + 2*nn_hls 
     
    8786         kpj = nn_jsize 
    8887      ENDIF 
    89   
     88      ! 
    9089      kpk = nn_ksize 
    9190      kperio = nn_perio 
    92  
    9391      !                             ! control print 
    94       ii = 1 
    95       WRITE(ldtxt(ii),*) '   '                                                                          ;   ii = ii + 1 
    96       WRITE(ldtxt(ii),*) 'usr_def_nam  : read the user defined namelist (namusr_def) in namelist_cfg'   ;   ii = ii + 1 
    97       WRITE(ldtxt(ii),*) '~~~~~~~~~~~ '                                                                 ;   ii = ii + 1 
    98       WRITE(ldtxt(ii),*) '   Namelist namusr_def : BENCH test case'                                     ;   ii = ii + 1 
    99       IF( nn_isize > 0 ) THEN 
    100          WRITE(ldtxt(ii),*) '      global domain size-x            nn_isize = ',  nn_isize              ;   ii = ii + 1 
    101       ELSE 
    102          WRITE(ldtxt(ii),*) '                                          jpni = ', jpni                   ;   ii = ii + 1 
    103          WRITE(ldtxt(ii),*) '       local domain size-x           -nn_isize = ', -nn_isize              ;   ii = ii + 1 
    104          WRITE(ldtxt(ii),*) '      global domain size-x                 kpi = ', kpi                    ;   ii = ii + 1 
     92      IF(lwp) THEN 
     93         WRITE(numout,*) '   ' 
     94         WRITE(numout,*) 'usr_def_nam  : read the user defined namelist (namusr_def) in namelist_cfg' 
     95         WRITE(numout,*) '~~~~~~~~~~~ ' 
     96         WRITE(numout,*) '   Namelist namusr_def : BENCH test case' 
     97         IF( nn_isize > 0 ) THEN 
     98            WRITE(numout,*) '      global domain size-x            nn_isize = ',  nn_isize 
     99         ELSE 
     100            WRITE(numout,*) '                                          jpni = ', jpni 
     101            WRITE(numout,*) '       local domain size-x           -nn_isize = ', -nn_isize 
     102            WRITE(numout,*) '      global domain size-x                 kpi = ', kpi 
     103         ENDIF 
     104         IF( nn_jsize > 0 ) THEN 
     105            WRITE(numout,*) '      global domain size-y            nn_jsize = ', nn_jsize 
     106         ELSE 
     107            WRITE(numout,*) '                                          jpnj = ', jpnj 
     108            WRITE(numout,*) '       local domain size-y           -nn_jsize = ', -nn_jsize 
     109            WRITE(numout,*) '      global domain size-y                 kpj = ', kpj 
     110         ENDIF 
     111         WRITE(numout,*) '      global domain size-z            nn_ksize = ', nn_ksize 
     112         WRITE(numout,*) '      LBC of the global domain          kperio = ', kperio 
    105113      ENDIF 
    106       IF( nn_jsize > 0 ) THEN 
    107          WRITE(ldtxt(ii),*) '      global domain size-y            nn_jsize = ', nn_jsize               ;   ii = ii + 1 
    108       ELSE 
    109          WRITE(ldtxt(ii),*) '                                          jpnj = ', jpnj                   ;   ii = ii + 1 
    110          WRITE(ldtxt(ii),*) '       local domain size-y           -nn_jsize = ', -nn_jsize              ;   ii = ii + 1 
    111          WRITE(ldtxt(ii),*) '      global domain size-y                 kpj = ', kpj                    ;   ii = ii + 1 
    112       ENDIF 
    113       WRITE(ldtxt(ii),*) '      global domain size-z            nn_ksize = ', nn_ksize                  ;   ii = ii + 1 
    114       WRITE(ldtxt(ii),*) '      LBC of the global domain          kperio = ', kperio                    ;   ii = ii + 1 
    115114      ! 
    116115   END SUBROUTINE usr_def_nam 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/tests/BENCH/MY_SRC/usrdef_zgr.F90

    r10170 r11586  
    194194      z2d(:,:) = REAL( jpkm1 , wp )          ! flat bottom 
    195195      ! 
     196      IF( jperio == 3 .OR. jperio ==4 ) THEN   ! add a small island in the upper corners to avoid model instabilities... 
     197         z2d(mi0(       1):mi1(     3),mj0(jpjglo-2):mj1(jpjglo)) = 0. 
     198         z2d(mi0(jpiglo-2):mi1(jpiglo),mj0(jpjglo-2):mj1(jpjglo)) = 0. 
     199      ENDIF 
     200      ! 
    196201      CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. )           ! set surrounding land to zero (here jperio=0 ==>> closed) 
    197202      ! 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/tests/BENCH/MY_SRC/zdfiwm.F90

    r10420 r11586  
    406406      REWIND( numnam_ref )              ! Namelist namzdf_iwm in reference namelist : Wave-driven mixing 
    407407      READ  ( numnam_ref, namzdf_iwm, IOSTAT = ios, ERR = 901) 
    408 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namzdf_iwm in reference namelist', lwp ) 
     408901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namzdf_iwm in reference namelist' ) 
    409409      ! 
    410410      REWIND( numnam_cfg )              ! Namelist namzdf_iwm in configuration namelist : Wave-driven mixing 
    411411      READ  ( numnam_cfg, namzdf_iwm, IOSTAT = ios, ERR = 902 ) 
    412 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namzdf_iwm in configuration namelist', lwp ) 
     412902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namzdf_iwm in configuration namelist' ) 
    413413      IF(lwm) WRITE ( numond, namzdf_iwm ) 
    414414      ! 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/tests/CANAL/EXPREF/context_nemo.xml

    r9930 r11586  
    66<context id="nemo"> 
    77<!-- $id$ --> 
     8    <variable_definition> 
     9    <!-- Year of time origin for NetCDF files; defaults to 1800 --> 
     10       <variable id="ref_year" type="int"   > 1800 </variable> 
     11       <variable id="rau0"     type="float" > 1026.0 </variable> 
     12       <variable id="cpocean"  type="float" > 3991.86795711963 </variable> 
     13       <variable id="convSpsu" type="float" > 0.99530670233846  </variable> 
     14       <variable id="rhoic"    type="float" > 917.0 </variable> 
     15       <variable id="rhosn"    type="float" > 330.0 </variable> 
     16       <variable id="missval"  type="float" > 1.e20 </variable> 
     17    </variable_definition> 
    818<!-- Fields definition --> 
    919    <field_definition src="./field_def_nemo-oce.xml"/>   <!--  NEMO ocean dynamics                     --> 
     
    1828     
    1929    <axis_definition> 
    20       <axis id="deptht" long_name="Vertical T levels" unit="m" positive="down" /> 
    21       <axis id="depthu" long_name="Vertical U levels" unit="m" positive="down" /> 
    22       <axis id="depthv" long_name="Vertical V levels" unit="m" positive="down" /> 
    23       <axis id="depthw" long_name="Vertical W levels" unit="m" positive="down" /> 
    24       <axis id="nfloat" long_name="Float number"      unit="-"                 /> 
    25       <axis id="icbcla"  long_name="Iceberg class"      unit="1"               /> 
    26       <axis id="ncatice" long_name="Ice category"       unit="1"               /> 
    27       <axis id="iax_20C" long_name="20 degC isotherm"   unit="degC"            /> 
    28       <axis id="iax_28C" long_name="28 degC isotherm"   unit="degC"            /> 
     30      <axis id="deptht"  long_name="Vertical T levels" unit="m" positive="down" /> 
     31      <axis id="depthu"  long_name="Vertical U levels" unit="m" positive="down" /> 
     32      <axis id="depthv"  long_name="Vertical V levels" unit="m" positive="down" /> 
     33      <axis id="depthw"  long_name="Vertical W levels" unit="m" positive="down" /> 
     34      <axis id="nfloat"  long_name="Float number"      unit="-"                 /> 
     35      <axis id="icbcla"  long_name="Iceberg class"     unit="1"                 /> 
     36      <axis id="ncatice" long_name="Ice category"      unit="1"                 /> 
     37      <axis id="iax_20C" long_name="20 degC isotherm"  unit="degC"              /> 
     38      <axis id="iax_28C" long_name="28 degC isotherm"  unit="degC"              /> 
    2939    </axis_definition> 
    3040  
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/tests/CANAL/EXPREF/namelist_cfg

    r10075 r11586  
    275275!!   namdiu       Cool skin and warm layer models                       (default: OFF) 
    276276!!   namdiu       Cool skin and warm layer models                       (default: OFF) 
    277 !!   namflo       float parameters                                      ("key_float") 
    278 !!   nam_diaharm  Harmonic analysis of tidal constituents               ("key_diaharm") 
    279 !!   namdct       transports through some sections                      ("key_diadct") 
     277!!   namflo       float parameters                                      (default: OFF) 
     278!!   nam_diaharm  Harmonic analysis of tidal constituents               (default: OFF) 
     279!!   nam_diadct   transports through some sections                      (default: OFF) 
    280280!!   nam_diatmb   Top Middle Bottom Output                              (default: OFF) 
    281281!!   nam_dia25h   25h Mean Output                                       (default: OFF) 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/tests/CANAL/MY_SRC/diawri.F90

    r10425 r11586  
    531531      !!      define all the NETCDF files and fields 
    532532      !!      At each time step call histdef to compute the mean if ncessary 
    533       !!      Each nwrite time step, output the instantaneous or mean fields 
     533      !!      Each nn_write time step, output the instantaneous or mean fields 
    534534      !!---------------------------------------------------------------------- 
    535535      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
     
    547547      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace 
    548548      !!---------------------------------------------------------------------- 
    549       !  
    550       IF( ln_timing )   CALL timing_start('dia_wri') 
    551549      ! 
    552550      IF( ninist == 1 ) THEN     !==  Output the initial state and forcings  ==! 
     
    555553      ENDIF 
    556554      ! 
     555      IF( nn_write == -1 )   RETURN   ! we will never do any output 
     556      !  
     557      IF( ln_timing )   CALL timing_start('dia_wri') 
     558      ! 
    557559      ! 0. Initialisation 
    558560      ! ----------------- 
     
    564566      clop = "x"         ! no use of the mask value (require less cpu time and otherwise the model crashes) 
    565567#if defined key_diainstant 
    566       zsto = nwrite * rdt 
     568      zsto = nn_write * rdt 
    567569      clop = "inst("//TRIM(clop)//")" 
    568570#else 
     
    570572      clop = "ave("//TRIM(clop)//")" 
    571573#endif 
    572       zout = nwrite * rdt 
     574      zout = nn_write * rdt 
    573575      zmax = ( nitend - nit000 + 1 ) * rdt 
    574576 
     
    601603         ! WRITE root name in date.file for use by postpro 
    602604         IF(lwp) THEN 
    603             CALL dia_nam( clhstnam, nwrite,' ' ) 
     605            CALL dia_nam( clhstnam, nn_write,' ' ) 
    604606            CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 
    605607            WRITE(inum,*) clhstnam 
     
    609611         ! Define the T grid FILE ( nid_T ) 
    610612 
    611          CALL dia_nam( clhstnam, nwrite, 'grid_T' ) 
     613         CALL dia_nam( clhstnam, nn_write, 'grid_T' ) 
    612614         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename 
    613615         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit 
     
    645647         ! Define the U grid FILE ( nid_U ) 
    646648 
    647          CALL dia_nam( clhstnam, nwrite, 'grid_U' ) 
     649         CALL dia_nam( clhstnam, nn_write, 'grid_U' ) 
    648650         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename 
    649651         CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu 
     
    658660         ! Define the V grid FILE ( nid_V ) 
    659661 
    660          CALL dia_nam( clhstnam, nwrite, 'grid_V' )                   ! filename 
     662         CALL dia_nam( clhstnam, nn_write, 'grid_V' )                   ! filename 
    661663         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam 
    662664         CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv 
     
    671673         ! Define the W grid FILE ( nid_W ) 
    672674 
    673          CALL dia_nam( clhstnam, nwrite, 'grid_W' )                   ! filename 
     675         CALL dia_nam( clhstnam, nn_write, 'grid_W' )                   ! filename 
    674676         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam 
    675677         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit 
     
    762764         ENDIF 
    763765 
    764          IF( .NOT. ln_cpl ) THEN 
     766         IF( ln_ssr ) THEN 
    765767            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp 
    766768               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
     
    770772               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    771773         ENDIF 
    772  
    773          IF( ln_cpl .AND. nn_ice <= 1 ) THEN 
    774             CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp 
    775                &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    776             CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp 
    777                &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    778             CALL histdef( nid_T, "sosafldp", "Surface salt flux: Damping"         , "Kg/m2/s",   &  ! erp * sn 
    779                &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    780          ENDIF 
    781           
     774        
    782775         clmx ="l_max(only(x))"    ! max index on a period 
    783776!         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX  
     
    855848      ! donne le nombre d'elements, et ndex la liste des indices a sortir 
    856849 
    857       IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN  
     850      IF( lwp .AND. MOD( itmod, nn_write ) == 0 ) THEN  
    858851         WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step' 
    859852         WRITE(numout,*) '~~~~~~ ' 
     
    919912      ENDIF 
    920913 
    921       IF( .NOT. ln_cpl ) THEN 
     914      IF( ln_ssr ) THEN 
    922915         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping 
    923916         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping 
    924          IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1) 
    925          CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping 
    926       ENDIF 
    927       IF( ln_cpl .AND. nn_ice <= 1 ) THEN 
    928          CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping 
    929          CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping 
    930          IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1) 
     917         zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1) 
    931918         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping 
    932919      ENDIF 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/tests/CANAL/MY_SRC/domvvl.F90

    r10425 r11586  
    994994      REWIND( numnam_ref )              ! Namelist nam_vvl in reference namelist :  
    995995      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901) 
    996 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp ) 
     996901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_vvl in reference namelist' ) 
    997997      REWIND( numnam_cfg )              ! Namelist nam_vvl in configuration namelist : Parameters of the run 
    998998      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 ) 
    999 902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp ) 
     999902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist' ) 
    10001000      IF(lwm) WRITE ( numond, nam_vvl ) 
    10011001      ! 
  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/tests/CANAL/MY_SRC/usrdef_nam.F90

    r10074 r11586  
    5858CONTAINS 
    5959 
    60    SUBROUTINE usr_def_nam( ldtxt, ldnam, cd_cfg, kk_cfg, kpi, kpj, kpk, kperio ) 
     60   SUBROUTINE usr_def_nam( cd_cfg, kk_cfg, kpi, kpj, kpk, kperio ) 
    6161      !!---------------------------------------------------------------------- 
    6262      !!                     ***  ROUTINE dom_nam  *** 
     
    7070      !! ** input   : - namusr_def namelist found in namelist_cfg 
    7171      !!---------------------------------------------------------------------- 
    72       CHARACTER(len=*), DIMENSION(:), INTENT(out) ::   ldtxt, ldnam    ! stored print information 
    7372      CHARACTER(len=*)              , INTENT(out) ::   cd_cfg          ! configuration name 
    7473      INTEGER                       , INTENT(out) ::   kk_cfg          ! configuration resolution 
     
    7675      INTEGER                       , INTENT(out) ::   kperio          ! lateral global domain b.c.  
    7776      ! 
    78       INTEGER ::   ios, ii      ! Local integer 
    79       REAL(wp)::   zh ! Local scalars 
     77      INTEGER ::   ios      ! Local integer 
     78      REAL(wp)::   zh       ! Local scalars 
    8079      !! 
    8180      NAMELIST/namusr_def/  rn_domszx, rn_domszy, rn_domszz, rn_dx, rn_dy, rn_dz, rn_0xratio, rn_0yratio   & 
     
    8584      !!---------------------------------------------------------------------- 
    8685      ! 
    87       ii = 1 
    88       ! 
    8986      REWIND( numnam_cfg )          ! Namelist namusr_def (exist in namelist_cfg only) 
    9087      READ  ( numnam_cfg, namusr_def, IOSTAT = ios, ERR = 902 ) 
    91 902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namusr_def in configuration namelist', .TRUE. ) 
     88902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namusr_def in configuration namelist', cdtxt ) 
     89      ! 
     90      IF(lwm)   WRITE( numond, namusr_def ) 
    9291      ! 
    9392#if defined key_agrif  
     
    103102#endif 
    104103      ! 
    105       WRITE( ldnam(:), namusr_def ) 
     104      IF(lwm)   WRITE( numond, namusr_def ) 
    106105      ! 
    107106      cd_cfg = 'EW_CANAL'             ! name & resolution (not used) 
     
    120119      ! 
    121120      zh  = (kpk-1)*rn_dz 
    122       !                             ! control print 
    123       WRITE(ldtxt(ii),*) '   '                                                                          ;   ii = ii + 1 
    124       WRITE(ldtxt(ii),*) 'usr_def_nam  : read the user defined namelist (namusr_def) in namelist_cfg'   ;   ii = ii + 1 
    125       WRITE(ldtxt(ii),*) '~~~~~~~~~~~ '                                                                 ;   ii = ii + 1 
    126       WRITE(ldtxt(ii),*) '   Namelist namusr_def : EW_CANAL test case'                                  ;   ii = ii + 1 
    127       WRITE(ldtxt(ii),*) '      horizontal domain size-x          rn_domszx  = ', rn_domszx, ' km'      ;   ii = ii + 1 
    128       WRITE(ldtxt(ii),*) '      horizontal domain size-y          rn_domszy  = ', rn_domszy, ' km'      ;   ii = ii + 1 
    129       WRITE(ldtxt(ii),*) '      vertical   domain size-z          rn_domszz  = ', rn_domszz, '  m'      ;   ii = ii + 1 
    130       WRITE(ldtxt(ii),*) '      horizontal x-resolution           rn_dx      = ',     rn_dx, ' km'      ;   ii = ii + 1 
    131       WRITE(ldtxt(ii),*) '      horizontal y-resolution           rn_dy      = ',     rn_dy, ' km'      ;   ii = ii + 1 
    132       WRITE(ldtxt(ii),*) '      vertical resolution               rn_dz      = ',     rn_dz, '  m'      ;   ii = ii + 1 
    133       WRITE(ldtxt(ii),*) '      x-domain ratio of the 0           rn_0xratio = ', rn_0xratio            ;   ii = ii + 1 
    134       WRITE(ldtxt(ii),*) '      y-domain ratio of the 0           rn_0yratio = ', rn_0yratio            ;   ii = ii + 1 
    135       WRITE(ldtxt(ii),*) '          H [m] : ', zh                                                       ;   ii = ii + 1 
    136       WRITE(ldtxt(ii),*) '      F computation                     nn_fcase   = ',   nn_fcase            ;   ii = ii + 1 
    137       WRITE(ldtxt(ii),*) '      Reference latitude                rn_ppgphi0 = ', rn_ppgphi0            ;   ii = ii + 1 
    138       WRITE(ldtxt(ii),*) '      10m wind speed                    rn_u10     = ',     rn_u10, ' m/s'    ;   ii = ii + 1 
    139       WRITE(ldtxt(ii),*) '         wind latitudinal extension     rn_windszy = ', rn_windszy, ' km'     ;   ii = ii + 1 
    140       WRITE(ldtxt(ii),*) '         wind longitudinal extension    rn_windszx = ', rn_windszx, ' km'     ;   ii = ii + 1 
    141       WRITE(ldtxt(ii),*) '         Uoce multiplicative factor     rn_uofac   = ',   rn_uofac            ;   ii = ii + 1 
    142       WRITE(ldtxt(ii),*) '      initial Canal max current        rn_vtxmax  = ',  rn_vtxmax, ' m/s'    ;   ii = ii + 1 
    143       WRITE(ldtxt(ii),*) '      initial zonal current             rn_uzonal  = ',  rn_uzonal, ' m/s'    ;   ii = ii + 1 
    144       WRITE(ldtxt(ii),*) '         Jet latitudinal extension      rn_ujetszy = ', rn_ujetszy, ' km'     ;   ii = ii + 1 
    145       WRITE(ldtxt(ii),*) '    &nbs