New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 11586 for NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src – NEMO

Ignore:
Timestamp:
2019-09-20T17:28:02+02:00 (5 years 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/src
Files:
37 edited

Legend:

Unmodified
Added
Removed

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