Ignore:
Timestamp:
2017-09-12T20:46:13+02:00 (3 years ago)
Author:
clem
Message:

changes in style - part6 - one more round

Location:
branches/2017/dev_r8183_ICEMODEL/NEMOGCM
Files:
17 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/ORCA2_SAS_LIM3/EXP00/namelist_ice_cfg

    r8321 r8517  
    11!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    2 !! LIM3 configuration namelist: Overwrites SHARED/namelist_ice_lim3_ref 
    3 !!              1 - Generic parameters                 (namicerun) 
    4 !!              2 - Diagnostics                        (namicediag) 
    5 !!              3 - Ice initialization                 (namiceini) 
    6 !!              4 - Ice discretization                 (namiceitd) 
    7 !!              5 - Ice dynamics and transport         (namicedyn) 
    8 !!              6 - Ice thermodynamics                 (namicethd) 
    9 !!              7 - Ice salinity                       (namicesal) 
    10 !!              8 - Ice mechanical redistribution      (namiceitdme) 
    11 !!              9 - Ice/snow albedos                   (namicealb) 
    12 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     2!! ESIM configuration namelist: Overwrites SHARED/namelist_ice_lim3_ref 
     3!!              1 - Generic parameters                 (namice_run) 
     4!!              2 - Ice thickness discretization       (namice_itd) 
     5!!              3 - Ice dynamics                       (namice_dyn) 
     6!!              4 - Ice ridging/rafting                (namice_rdgrft) 
     7!!              5 - Ice rheology                       (namice_rhg) 
     8!!              6 - Ice advection                      (namice_adv) 
     9!!              7 - Ice thermodynamics                 (namice_thd) 
     10!!              8 - Ice salinity                       (namice_sal) 
     11!!              9 - Ice melt ponds                     (namice_mp) 
     12!!             10 - Ice initialization                 (namice_ini) 
     13!!             11 - Ice/snow albedos                   (namice_alb) 
     14!!             12 - Ice diagnostics                    (namice_dia) 
     15!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     16! 
    1317!------------------------------------------------------------------------------ 
    14 &namicerun     !   Generic parameters 
     18&namice_run     !   Generic parameters 
    1519!------------------------------------------------------------------------------ 
    1620/ 
    1721!------------------------------------------------------------------------------ 
    18 &namicediag    !   Diagnostics 
     22&namice_itd     !   Ice discretization 
    1923!------------------------------------------------------------------------------ 
    2024/ 
    2125!------------------------------------------------------------------------------ 
    22 &namiceini     !   Ice initialization 
     26&namice_dyn     !   Ice dynamics 
    2327!------------------------------------------------------------------------------ 
    2428/ 
    2529!------------------------------------------------------------------------------ 
    26 &namiceitd     !   Ice discretization 
     30&namice_rdgrft  !   Ice ridging/rafting 
    2731!------------------------------------------------------------------------------ 
    2832/ 
    2933!------------------------------------------------------------------------------ 
    30 &namicedyn     !   Ice dynamics and transport 
     34&namice_rhg     !   Ice rheology 
    3135!------------------------------------------------------------------------------ 
    3236/ 
    3337!------------------------------------------------------------------------------ 
    34 &namicethd     !   Ice thermodynamics 
     38&namice_adv     !   Ice advection 
    3539!------------------------------------------------------------------------------ 
    3640/ 
    3741!------------------------------------------------------------------------------ 
    38 &namicesal     !   Ice salinity 
     42&namice_thd     !   Ice thermodynamics 
    3943!------------------------------------------------------------------------------ 
    4044/ 
    4145!------------------------------------------------------------------------------ 
    42 &namiceitdme   !   Ice mechanical redistribution (ridging and rafting) 
     46&namice_sal     !   Ice salinity 
    4347!------------------------------------------------------------------------------ 
    4448/ 
    45 !----------------------------------------------------------------------- 
    46 &namicealb     !   albedo parameters 
    47 !----------------------------------------------------------------------- 
     49!------------------------------------------------------------------------------ 
     50&namicemp      !   Melt ponds 
     51!------------------------------------------------------------------------------ 
    4852/ 
     53!------------------------------------------------------------------------------ 
     54&namice_ini     !   Ice initialization 
     55!------------------------------------------------------------------------------ 
     56/ 
     57!------------------------------------------------------------------------------ 
     58&namice_alb     !   albedo parameters 
     59!------------------------------------------------------------------------------ 
     60/ 
     61!------------------------------------------------------------------------------ 
     62&namice_dia     !   Diagnostics 
     63!------------------------------------------------------------------------------ 
     64/ 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref

    r8516 r8517  
    1818&namice_run     !   Generic parameters 
    1919!------------------------------------------------------------------------------ 
    20    jpl              =    5            !  number of ice  categories 
    21    nlay_i           =    2            !  number of ice  layers 
    22    nlay_s           =    1            !  number of snow layers (only 1 is working) 
    23    nn_monocat       =    0            !  virtual ITD mono-category parameterizations (1-4 => jpl = 1 only) or not (0) 
     20   jpl              =   5             !  number of ice  categories 
     21   nlay_i           =   2             !  number of ice  layers 
     22   nlay_s           =   1             !  number of snow layers (only 1 is working) 
     23   nn_monocat       =   0             !  virtual ITD mono-category parameterizations (1-4 => jpl = 1 only) or not (0) 
    2424                                      !     2: simple piling instead of ridging    --- temporary option 
    2525                                      !     3: activate G(he) only                 --- temporary option 
    2626                                      !     4: activate extra lateral melting only --- temporary option 
     27   ln_icedyn        = .true.          !  ice dynamics (T) or not (F) 
     28   ln_icethd        = .true.          !  ice thermo   (T) or not (F) 
    2729   rn_amax_n        =   0.997         !  maximum tolerated ice concentration NH 
    2830   rn_amax_s        =   0.997         !  maximum tolerated ice concentration SH 
     
    4143&namice_dyn     !   Ice dynamics 
    4244!------------------------------------------------------------------------------ 
    43    ln_icedyn        = .true.          !  ice dynamics (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 
    44       nn_icedyn     =   2             !     switch for ice dynamics    
    45                                       !        2: total 
    46                                       !        1: advection only (no diffusion, no ridging/rafting) 
    47                                       !        0: advection only (as 1 but with prescribed velocity, bypass rheology) 
    48          rn_uice    =   0.00001       !           prescribed ice u-velocity 
    49          rn_vice    =  -0.00001       !           prescribed ice v-velocity 
     45   ln_dynFULL       = .true.          !  dyn.: full ice dynamics               (rheology + advection + ridging/rafting + correction) 
     46   ln_dynRHGADV     = .false.         !  dyn.: no ridge/raft & no corrections  (rheology + advection) 
     47   ln_dynADV        = .false.         !  dyn.: only advection w prescribed vel.(rn_uvice + advection) 
     48      rn_uice       =   0.00001       !        prescribed ice u-velocity 
     49      rn_vice       =   0.            !        prescribed ice v-velocity 
    5050   rn_ishlat        =   2.            !  free slip (0) ; partial slip (0-2) ; no slip (2) ; strong slip (>2) 
    5151   rn_cio           =   5.0e-03       !  ice-ocean drag coefficient (-) 
     
    5353      rn_gamma      =   0.15          !     fraction of ocean depth that ice must reach to initiate landfast 
    5454                                      !        recommended range: [0.1 ; 0.25] 
    55       rn_icebfr     =  10.            !     maximum bottom stress per unit area of contact (N/m2)                  
     55      rn_icebfr     =  10.            !     maximum bottom stress per unit area of contact [N/m2]                  
    5656                                      !        a very large value ensures ice velocity=0 even with a small contact area 
    5757                                      !        recommended range: ?? (should be greater than atm-ice stress => >0.1 N/m2) 
    58       rn_lfrelax    =   1.e-5         !     relaxation time scale to reach static friction (s-1)                  
     58      rn_lfrelax    =   1.e-5         !     relaxation time scale to reach static friction [s-1] 
    5959/ 
    6060!------------------------------------------------------------------------------ 
     
    6363          ! -- ice_rdgrft_strength -- ! 
    6464   ln_str_H79       = .true.          !  ice strength param.: Hibler_79   => P = pstar*<h>*exp(-c_rhg*A)                       
    65       rn_pstar      =   2.0e+04       !     ice strength thickness parameter (N/m2)  
     65      rn_pstar      =   2.0e+04       !     ice strength thickness parameter [N/m2] 
    6666      rn_crhg       =   20.0          !     ice strength conc. parameter (-) 
    6767   ln_str_R75       = .false.         !  ice strength param.: Rothrock_75 => P = Cf*coeff*integral(wr.h^2)     
     
    7474   ln_partf_exp     = .true.          !  Exponential ridging participation function (Lipscomb, 2007) 
    7575      rn_astar      =   0.03          !     exponential measure of ridging ice fraction [set to 0.05 if hstar=100] 
    76    ln_ridging       = .true.          !  ridging activated (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 
    77       rn_hstar      =  25.0           !     determines the maximum thickness of ridged ice (m) (Hibler, 1980) 
     76   ln_ridging       = .true.          !  ridging activated (T) or not (F) 
     77      rn_hstar      =  25.0           !     determines the maximum thickness of ridged ice [m] (Hibler, 1980) 
    7878      rn_porordg    =   0.3           !     porosity of newly ridged ice (Lepparanta et al., 1995) 
    7979      rn_fsnwrdg    =   0.5           !     snow volume fraction that survives in ridging 
    8080      rn_fpndrdg    =   1.0           !     pond fraction that survives in ridging (small a priori) 
    81    ln_rafting       = .true.          !  rafting activated (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 
    82       rn_hraft      =   0.75          !     threshold thickness for rafting (m) 
     81   ln_rafting       = .true.          !  rafting activated (T) or not (F) 
     82      rn_hraft      =   0.75          !     threshold thickness for rafting [m] 
    8383      rn_craft      =   5.0           !     squeezing coefficient used in the rafting function 
    8484      rn_fsnwrft    =   0.5           !     snow volume fraction that survives in rafting 
     
    8989!------------------------------------------------------------------------------ 
    9090   ln_rhg_EVP       = .true.          !  EVP rheology 
    91       rn_creepl     =   1.0e-12       !     creep limit (s-1) 
     91      rn_creepl     =   1.0e-12       !     creep limit [1/s] 
    9292      rn_ecc        =   2.0           !     eccentricity of the elliptical yield curve           
    9393      nn_nevp       = 120             !     number of EVP subcycles                              
     
    105105&namice_thd     !   Ice thermodynamics 
    106106!------------------------------------------------------------------------------ 
    107    ln_icethd        = .true.            !  ice thermo   (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 
    108107                     ! -- icethd_dif -- ! 
    109    rn_kappa_i       =   1.0             !  radiation attenuation coefficient in sea ice (m-1) 
     108   rn_kappa_i       =   1.0             !  radiation attenuation coefficient in sea ice [1/m] 
    110109   ln_cndi_U64      = .false.           !  sea ice thermal conductivity: k = k0 + beta.S/T            (Untersteiner, 1964) 
    111110   ln_cndi_P07      = .true.            !  sea ice thermal conductivity: k = k0 + beta1.S/T - beta2.T (Pringle et al., 2007) 
     
    114113                                        !     Obs: 0.1-0.5 (Lecomte et al, JAMES 2013) 
    115114                      ! -- icethd_dh -- ! 
    116    ln_icedH         = .true.            !  activate ice thickness change from growing/melting (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 
     115   ln_icedH         = .true.            !  activate ice thickness change from growing/melting (T) or not (F) 
    117116   rn_blow_s        =   0.66            !     mesure of snow blowing into the leads 
    118117                                        !     = 1 => no snow blowing, < 1 => some snow blowing 
    119118                      ! -- icethd_da -- ! 
    120    ln_icedA         = .true.            !  activate lateral melting param. (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 
     119   ln_icedA         = .true.            !  activate lateral melting param. (T) or not (F) 
    121120      rn_beta       =   1.0             !     coef. beta for lateral melting param. Recommended range=[0.8-1.2] 
    122121                                        !      => decrease = more melt and melt peaks toward higher concentration (A~0.5 for beta=1 ; A~0.8 for beta=0.2) 
     
    127126                                        !         10 vs 8m = -20% melting 
    128127                     ! -- icethd_lac -- ! 
    129    ln_icedO         = .true.            !  activate ice growth in open-water (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 
     128   ln_icedO         = .true.            !  activate ice growth in open-water (T) or not (F) 
    130129      rn_hinew      =   0.1             !     thickness for new ice formation in open water (m), must be larger than rn_hnewice 
    131130      ln_frazil     = .false.           !     Frazil ice parameterization (ice collection as a function of wind) 
     
    148147!------------------------------------------------------------------------------ 
    149148                    ! -- icethd_sal -- ! 
    150    ln_icedS       = .true.             !  activate gravity drainage and flushing (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 
     149   ln_icedS       = .true.             !  activate gravity drainage and flushing (T) or not (F) 
    151150   nn_icesal      =   2                !  ice salinity option 
    152151                                       !     1: constant ice salinity (S=rn_icesal) 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ref

    r8321 r8517  
    429429/ 
    430430!----------------------------------------------------------------------- 
    431 &namsbc_alb    !   albedo parameters 
    432 !----------------------------------------------------------------------- 
    433    nn_ice_alb   =    1   !  parameterization of ice/snow albedo 
    434                          !     0: Shine & Henderson-Sellers (JGR 1985), giving clear-sky albedo 
    435                          !     1: "home made" based on Brandt et al. (JClim 2005) and Grenfell & Perovich (JGR 2004), 
    436                          !        giving cloud-sky albedo 
    437    rn_alb_sdry  =  0.85  !  dry snow albedo         : 0.80 (nn_ice_alb = 0); 0.85 (nn_ice_alb = 1); obs 0.85-0.87 (cloud-sky) 
    438    rn_alb_smlt  =  0.75  !  melting snow albedo     : 0.65 ( '' )          ; 0.75 ( '' )          ; obs 0.72-0.82 ( '' ) 
    439    rn_alb_idry  =  0.60  !  dry ice albedo          : 0.72 ( '' )          ; 0.60 ( '' )          ; obs 0.54-0.65 ( '' ) 
    440    rn_alb_imlt  =  0.50  !  bare puddled ice albedo : 0.53 ( '' )          ; 0.50 ( '' )          ; obs 0.49-0.58 ( '' ) 
    441    rn_alb_dpnd  =  0.27  !  ponded ice albedo       : 0.25 ( '' )          ; 0.27 ( '' )          ; obs 0.10-0.30 ( '' ) 
    442 / 
    443 !----------------------------------------------------------------------- 
    444431&namsbc_wave   ! External fields from wave model                        (ln_wave=T) 
    445432!----------------------------------------------------------------------- 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r8516 r8517  
    157157   INTEGER           , PUBLIC ::   nlay_i          !: number of ice  layers  
    158158   INTEGER           , PUBLIC ::   nlay_s          !: number of snow layers  
     159   LOGICAL           , PUBLIC ::   ln_icedyn       !: flag for ice dynamics (T) or not (F) 
     160   LOGICAL           , PUBLIC ::   ln_icethd       !: flag for ice thermo   (T) or not (F) 
    159161   REAL(wp)          , PUBLIC ::   rn_amax_n       !: maximum ice concentration Northern hemisphere 
    160162   REAL(wp)          , PUBLIC ::   rn_amax_s       !: maximum ice concentration Southern hemisphere 
     
    168170    
    169171   !                                     !!** ice-dynamics namelist (namice_dyn) ** 
    170    LOGICAL , PUBLIC ::   ln_icedyn        !: flag for ice dynamics (T) or not (F) 
    171    INTEGER , PUBLIC ::   nn_icedyn        !: flag for ice dynamics 
    172    REAL(wp), PUBLIC ::   rn_uice          !: prescribed u-vel (case nn_icedyn=0) 
    173    REAL(wp), PUBLIC ::   rn_vice          !: prescribed v-vel (case nn_icedyn=0) 
     172   LOGICAL,  PUBLIC ::   ln_dynFULL       !:  dyn.: full ice dynamics               (rheology + advection + ridging/rafting + correction) 
     173   LOGICAL,  PUBLIC ::   ln_dynRHGADV     !:  dyn.: no ridge/raft & no corrections  (rheology + advection) 
     174   LOGICAL,  PUBLIC ::   ln_dynADV        !:  dyn.: only advection w prescribed vel.(rn_uvice + advection) 
    174175   REAL(wp), PUBLIC ::   rn_ishlat        !: lateral boundary condition for sea-ice 
    175176   REAL(wp), PUBLIC ::   rn_cio           !: drag coefficient for oceanic stress 
    176177   LOGICAL , PUBLIC ::   ln_landfast      !: landfast ice parameterization (T or F)  
    177    REAL(wp), PUBLIC ::   rn_gamma         !: fraction of ocean depth that ice must reach to initiate landfast ice 
    178    REAL(wp), PUBLIC ::   rn_icebfr        !: maximum bottom stress per unit area of contact (landfast ice)  
    179    REAL(wp), PUBLIC ::   rn_lfrelax       !: relaxation time scale (s-1) to reach static friction (landfast ice)  
     178   REAL(wp), PUBLIC ::   rn_gamma         !:    fraction of ocean depth that ice must reach to initiate landfast ice 
     179   REAL(wp), PUBLIC ::   rn_icebfr        !:    maximum bottom stress per unit area of contact (landfast ice)  
     180   REAL(wp), PUBLIC ::   rn_lfrelax       !:    relaxation time scale (s-1) to reach static friction (landfast ice)  
    180181   ! 
    181182   !                                     !!** ice-rdige/raft namelist (namice_rdgrft) ** 
     
    199200   ! 
    200201   !                                     !!** ice-thermodynamics namelist (namice_thd) ** 
    201    LOGICAL , PUBLIC ::   ln_icethd        !: flag for ice thermo (T) or not (F) 
    202202                                          ! -- icethd_dif -- ! 
    203203   REAL(wp), PUBLIC ::   rn_kappa_i       !: coef. for the extinction of radiation Grenfell et al. (2006) [1/m] 
     
    266266   REAL(wp), PUBLIC ::   r1_nlay_s        !: 1 / nlay_s  
    267267   REAL(wp), PUBLIC ::   rswitch          !: switch for the presence of ice (1) or not (0) 
     268   REAL(wp), PUBLIC ::   rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft   !: conservation diagnostics 
    268269   REAL(wp), PUBLIC, PARAMETER ::   epsi06 = 1.e-06_wp  !: small number  
    269270   REAL(wp), PUBLIC, PARAMETER ::   epsi10 = 1.e-10_wp  !: small number  
    270271   REAL(wp), PUBLIC, PARAMETER ::   epsi20 = 1.e-20_wp  !: small number  
    271    ! 
    272    LOGICAL , PUBLIC ::   l_piling         !: =T simple conservative piling, comparable with LIM2 
     272 
    273273 
    274274   !                                     !!** define arrays 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90

    r8516 r8517  
    3737   PUBLIC   ice_adv_init   ! called by icestp 
    3838 
     39   INTEGER ::              nice_adv   ! choice of the type of advection scheme 
     40   !                                        ! associated indices: 
     41   INTEGER, PARAMETER ::   np_advPRA = 1   ! Prather scheme 
     42   INTEGER, PARAMETER ::   np_advUMx = 2   ! Ultimate-Macho scheme 
     43 
    3944   !! * Substitution 
    4045#  include "vectopt_loop_substitute.h90" 
     
    6166      !!---------------------------------------------------------------------- 
    6267      INTEGER, INTENT(in) ::   kt   ! number of iteration 
    63       ! 
    64       INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices 
    65       INTEGER  ::   initad                  ! number of sub-timestep for the advection 
    66       REAL(wp) ::   zcfl , zusnit           !   -      - 
    67       CHARACTER(len=80) :: cltmp 
    68       ! 
    69       REAL(wp) ::    zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
    70       REAL(wp) ::    zdv 
    71       REAL(wp), DIMENSION(jpi,jpj)           ::   zatold, zeiold, zesold, zsmvold  
    72       REAL(wp), DIMENSION(jpi,jpj,jpl)       ::   zhimax, zviold, zvsold 
    7368      !!--------------------------------------------------------------------- 
    7469      IF( nn_timing == 1 )  CALL timing_start('iceadv') 
     
    8075      ENDIF 
    8176       
    82       CALL ice_var_agg( 1 ) ! integrated values + ato_i 
    83  
    8477      ! conservation test 
    85       IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceadv', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    86        
    87       ! store old values for diag 
    88       zviold (:,:,:) = v_i(:,:,:) 
    89       zvsold (:,:,:) = v_s(:,:,:) 
    90       zsmvold(:,:)   = SUM( smv_i(:,:,:), dim=3 ) 
    91       zeiold (:,:)   = et_i(:,:) 
    92       zesold (:,:)   = et_s(:,:) 
    93  
    94       ! Thickness correction init. 
    95       zatold(:,:) = at_i 
    96       WHERE( a_i(:,:,:) >= epsi20 ) 
    97          ht_i(:,:,:) = v_i(:,:,:) / a_i(:,:,:) 
    98          ht_s(:,:,:) = v_s(:,:,:) / a_i(:,:,:) 
    99       ELSEWHERE 
    100          ht_i(:,:,:) = 0._wp 
    101          ht_s(:,:,:) = 0._wp          
    102       END WHERE 
    103           
    104       ! Record max of the surrounding ice thicknesses for correction in case advection creates ice too thick 
    105       zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:) 
    106       DO jl = 1, jpl 
    107          DO jj = 2, jpjm1 
    108             DO ji = 2, jpim1 
    109 !!gm use of MAXVAL here is very probably less efficient than expending the 9 values 
    110                zhimax(ji,jj,jl) = MAX( epsi20, MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) + ht_s(ji-1:ji+1,jj-1:jj+1,jl) ) ) 
    111             END DO 
    112          END DO 
    113       END DO 
    114       CALL lbc_lnk( zhimax(:,:,:), 'T', 1. ) 
    115        
     78      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceadv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
     79                      
    11680      !---------- 
    11781      ! Advection 
    11882      !---------- 
    119       IF( ln_adv_UMx ) THEN                    !-- ULTIMATE-MACHO scheme 
     83      SELECT CASE( nice_adv ) 
     84      !                                !-----------------------! 
     85      CASE( np_advUMx )                ! ULTIMATE-MACHO scheme ! 
     86         !                             !-----------------------! 
    12087         CALL ice_adv_umx( kt, u_ice, v_ice,  & 
    12188            &              ato_i, v_i, v_s, smv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
    122           
    123       ELSEIF( ln_adv_Pra ) THEN                !-- PRATHER scheme 
     89      !                                !-----------------------! 
     90      CASE( np_advPRA )                ! PRATHER scheme        ! 
     91         !                             !-----------------------! 
    12492         CALL ice_adv_prather( kt, u_ice, v_ice,  & 
    12593            &                  ato_i, v_i, v_s, smv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
    126           
    127       ENDIF 
     94      END SELECT 
    12895 
    129       ! total ice fraction 
    130       at_i(:,:) = a_i(:,:,1) 
    131       DO jl = 2, jpl 
    132          at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 
    133       END DO 
    134        
    13596      !------------ 
    13697      ! diagnostics 
    13798      !------------ 
    138       DO jj = 1, jpj 
    139          DO ji = 1, jpi 
    140             diag_trp_ei (ji,jj) = ( SUM( e_i  (ji,jj,1:nlay_i,:) ) -  zeiold(ji,jj) ) * r1_rdtice 
    141             diag_trp_es (ji,jj) = ( SUM( e_s  (ji,jj,1:nlay_s,:) ) -  zesold(ji,jj) ) * r1_rdtice 
    142             diag_trp_smv(ji,jj) = ( SUM( smv_i(ji,jj,:)          ) - zsmvold(ji,jj) ) * r1_rdtice 
    143             diag_trp_vi (ji,jj) =   SUM(   v_i(ji,jj,:)            -  zviold(ji,jj,:) ) * r1_rdtice 
    144             diag_trp_vs (ji,jj) =   SUM(   v_s(ji,jj,:)            -  zvsold(ji,jj,:) ) * r1_rdtice 
    145          END DO 
    146       END DO 
     99      diag_trp_ei (:,:) = SUM(SUM( e_i  (:,:,1:nlay_i,:) - e_i_b  (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_rdtice 
     100      diag_trp_es (:,:) = SUM(SUM( e_s  (:,:,1:nlay_s,:) - e_s_b  (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_rdtice 
     101      diag_trp_smv(:,:) = SUM(     smv_i(:,:,:)          - smv_i_b(:,:,:)                  , dim=3 ) * r1_rdtice 
     102      diag_trp_vi (:,:) = SUM(     v_i  (:,:,:)          - v_i_b  (:,:,:)                  , dim=3 ) * r1_rdtice 
     103      diag_trp_vs (:,:) = SUM(     v_s  (:,:,:)          - v_s_b  (:,:,:)                  , dim=3 ) * r1_rdtice 
    147104      IF( iom_use('icetrp') )   CALL iom_put( "icetrp" , diag_trp_vi * rday  )         ! ice volume transport 
    148105      IF( iom_use('snwtrp') )   CALL iom_put( "snwtrp" , diag_trp_vs * rday  )         ! snw volume transport 
     
    150107      IF( iom_use('deitrp') )   CALL iom_put( "deitrp" , diag_trp_ei         )         ! advected ice enthalpy (W/m2) 
    151108      IF( iom_use('destrp') )   CALL iom_put( "destrp" , diag_trp_es         )         ! advected snw enthalpy (W/m2) 
    152        
    153       !-------------------------------------- 
    154       ! Thickness correction in case too high 
    155       !-------------------------------------- 
    156       IF( nn_icedyn == 2 ) THEN 
    157          ! 
    158          CALL ice_var_zapsmall                       !-- zap small areas 
    159          ! 
    160          DO jl = 1, jpl 
    161             DO jj = 1, jpj 
    162                DO ji = 1, jpi 
    163                   IF ( v_i(ji,jj,jl) > 0._wp ) THEN  !-- bound to zhimax 
    164                      ! 
    165                      ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl) 
    166                      ht_s  (ji,jj,jl) = v_s (ji,jj,jl) / a_i(ji,jj,jl) 
    167                      zdv = v_i(ji,jj,jl) + v_s(ji,jj,jl) - zviold(ji,jj,jl) - zvsold(ji,jj,jl)   
    168                      ! 
    169                      IF ( ( zdv >  0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) .AND. zatold(ji,jj) < 0.80 ) .OR. & 
    170                         & ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN 
    171                         a_i (ji,jj,jl) = ( v_i(ji,jj,jl) + v_s(ji,jj,jl) ) / zhimax(ji,jj,jl) 
    172                         ht_i(ji,jj,jl) =   v_i(ji,jj,jl) / a_i(ji,jj,jl) 
    173                      ENDIF 
    174                      ! 
    175                   ENDIF 
    176                END DO 
    177             END DO 
    178          END DO 
    179                    
    180          WHERE( ht_i(:,:,jpl) > hi_max(jpl) )        !-- bound ht_i to hi_max (99 m) 
    181             ht_i(:,:,jpl) = hi_max(jpl) 
    182             a_i (:,:,jpl) = v_i(:,:,jpl) / hi_max(jpl) 
    183          END WHERE 
    184  
    185          IF ( nn_pnd_scheme > 0 ) THEN               !-- correct pond fraction to avoid a_ip > a_i 
    186             WHERE( a_ip(:,:,:) > a_i(:,:,:) )   a_ip(:,:,:) = a_i(:,:,:) 
    187          ENDIF 
    188          ! 
    189       ENDIF 
    190           
    191       !------------------------------------------------------------ 
    192       ! Impose a_i < amax if no ridging/rafting or in mono-category 
    193       !------------------------------------------------------------ 
    194       IF( l_piling ) THEN                            !-- simple conservative piling, comparable with 1-cat models 
    195          at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 
    196          DO jl = 1, jpl 
    197             WHERE( at_i(:,:) > epsi20 ) 
    198                a_i(:,:,jl) = a_i(:,:,jl) * (  1._wp + MIN( rn_amax_2d(:,:) - at_i(:,:) , 0._wp ) / at_i(:,:)  ) 
    199             END WHERE 
    200          END DO 
    201       ENDIF 
    202        
    203       ! agglomerate variables  
    204       vt_i(:,:) = SUM( v_i(:,:,:), dim=3 ) 
    205       vt_s(:,:) = SUM( v_s(:,:,:), dim=3 ) 
    206       at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 
    207  
    208       ! MV MP 2016 (remove once we get rid of a_i_frac and ht_i) 
    209       IF ( nn_pnd_scheme > 0 ) THEN 
    210           at_ip(:,:) = SUM( a_ip(:,:,:), dim = 3 ) 
    211           vt_ip(:,:) = SUM( v_ip(:,:,:), dim = 3 ) 
    212       ENDIF 
    213       ! END MP 2016 
    214        
    215       ! open water = 1 if at_i=0  
    216       WHERE( at_i == 0._wp ) ato_i = 1._wp  
    217        
    218       ! conservation test 
    219       IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceadv', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    220          
    221       ! -------------- 
    222       ! control prints 
    223       ! -------------- 
    224       IF( ln_icectl )   CALL ice_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' ) 
    225       ! 
    226       IF( nn_timing == 1 )  CALL timing_stop('iceadv') 
     109                
     110      ! controls 
     111      IF( ln_icediachk   )   CALL ice_cons_hsm(1, 'iceadv', 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 dyn & trp - ')                       ! prints 
     113      IF( nn_timing == 1 )   CALL timing_stop ('iceadv')                                                             ! timing 
    227114      ! 
    228115   END SUBROUTINE ice_adv 
     
    241128      !! ** input   :   Namelist namice_adv 
    242129      !!------------------------------------------------------------------- 
    243       INTEGER ::   ios   ! Local integer output status for namelist read 
     130      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read 
    244131      !! 
    245132      NAMELIST/namice_adv/ ln_adv_Pra, ln_adv_UMx, nn_UMx 
     
    260147         WRITE(numout,*) '~~~~~~~~~~~~' 
    261148         WRITE(numout,*) '   Namelist namice_adv' 
    262          WRITE(numout,*) '      advection scheme for ice transport (limtrp)' 
    263          WRITE(numout,*) '         type of advection scheme (Prather)                     ln_adv_Pra = ', ln_adv_Pra  
    264          WRITE(numout,*) '         type of advection scheme (Ulimate-Macho)               ln_adv_UMx = ', ln_adv_UMx  
    265          WRITE(numout,*) '            order of the Ultimate-Macho scheme                      nn_UMx = ', nn_UMx 
     149         WRITE(numout,*) '      type of advection scheme (Prather)                     ln_adv_Pra = ', ln_adv_Pra  
     150         WRITE(numout,*) '      type of advection scheme (Ulimate-Macho)               ln_adv_UMx = ', ln_adv_UMx  
     151         WRITE(numout,*) '         order of the Ultimate-Macho scheme                      nn_UMx = ', nn_UMx 
    266152      ENDIF 
    267153      ! 
    268       IF ( ( ln_adv_Pra .AND. ln_adv_UMx ) .OR. ( .NOT.ln_adv_Pra .AND. .NOT.ln_adv_UMx ) ) THEN 
    269          CALL ctl_stop( 'ice_adv_init: choose one and only one ice advection scheme (ln_adv_Pra or ln_adv_UMx)' ) 
    270       ENDIF 
     154      !                             !== set the choice of ice advection ==! 
     155      ioptio = 0  
     156      IF( ln_adv_Pra ) THEN   ;   ioptio = ioptio + 1   ;   nice_adv = np_advPRA    ;   ENDIF 
     157      IF( ln_adv_UMx ) THEN   ;   ioptio = ioptio + 1   ;   nice_adv = np_advUMx    ;   ENDIF 
     158      IF( ioptio /= 1 )   CALL ctl_stop( 'ice_adv_init: choose one and only one ice advection scheme (ln_adv_Pra or ln_adv_UMx)' ) 
    271159      ! 
    272160   END SUBROUTINE ice_adv_init 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90

    r8514 r8517  
    321321901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_alb in reference namelist', lwp ) 
    322322 
    323       REWIND( numnam_ice_cfg )              ! Namelist namsbc_alb in configuration namelist : Albedo parameters 
     323      REWIND( numnam_ice_cfg )              ! Namelist namice_alb in configuration namelist : Albedo parameters 
    324324      READ  ( numnam_ice_cfg, namice_alb, IOSTAT = ios, ERR = 902 ) 
    325325902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_alb in configuration namelist', lwp ) 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icecor.F90

    r8514 r8517  
    5353      ! 
    5454      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    55       REAL(wp) ::   zsal, zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b, zzc 
     55      REAL(wp) ::   zsal, zzc 
    5656      REAL(wp), DIMENSION(jpi,jpj) ::   zafx   ! concentration trends diag 
    5757      !!---------------------------------------------------------------------- 
    58       IF( nn_timing == 1 )   CALL timing_start('icecor') 
     58      ! controls 
     59      IF( nn_timing == 1 )   CALL timing_start('icecor')                                                             ! timing 
     60      IF( ln_icediachk   )   CALL ice_cons_hsm(0, 'icecor', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    5961      ! 
    6062      IF( kt == nit000 .AND. lwp .AND. kn == 2 ) THEN 
     
    6365         WRITE(numout,*) '~~~~~~~' 
    6466      ENDIF 
    65       !                             !--- conservation test 
    66       IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icecor', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    67       ! 
    6867      ! 
    6968      IF( kn == 2 ) THEN 
     
    175174      END SELECT 
    176175      ! 
    177       !                                !--- conservation test 
    178       IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icecor', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    179       ! 
    180       !                                !--- control prints 
    181       IF( ln_ctl )                    CALL ice_prt3D( 'icecor' ) 
    182       IF( ln_icectl .AND. kn == 2 )   CALL ice_prt( kt, iiceprt, jiceprt, 2, ' - Final state - ' ) 
    183       ! 
    184       IF( nn_timing == 1 )   CALL timing_stop('icecor') 
     176      ! controls 
     177      IF( ln_icediachk   )   CALL ice_cons_hsm(1, 'icecor', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
     178      IF( ln_ctl         )   CALL ice_prt3D   ('icecor')                                                             ! prints 
     179      IF( ln_icectl .AND. kn == 2 )   CALL ice_prt( kt, iiceprt, jiceprt, 2, ' - Final state - ' )                   ! prints 
     180      IF( nn_timing == 1 )   CALL timing_stop ('icecor')                                                             ! timing 
    185181      ! 
    186182   END SUBROUTINE ice_cor 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90

    r8514 r8517  
    4747CONTAINS 
    4848 
    49    SUBROUTINE ice_cons_hsm( icount, cd_routine, zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b ) 
     49   SUBROUTINE ice_cons_hsm( icount, cd_routine, pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft ) 
    5050      !!---------------------------------------------------------------------- 
    5151      !!                       ***  ROUTINE ice_cons_hsm *** 
     
    5656      !! ** Method  : This is an online diagnostics which can be activated with ln_icediachk=true 
    5757      !!              It prints in ocean.output if there is a violation of conservation at each time-step 
    58       !!              The thresholds (zv_sill, zs_sill, zh_sill) which determine violations are set to 
     58      !!              The thresholds (zv_sill, zs_sill, zt_sill) which determine violations are set to 
    5959      !!              a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years. 
    6060      !!              For salt and heat thresholds, ice is considered to have a salinity of 10  
     
    6363      INTEGER         , INTENT(in)    ::   icount        ! called at: =0 the begining of the routine, =1  the end 
    6464      CHARACTER(len=*), INTENT(in)    ::   cd_routine    ! name of the routine 
    65       REAL(wp)        , INTENT(inout) ::   zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b   ! ???? 
    66       !! 
    67       REAL(wp)                        :: zvi,   zsmv,   zei,   zfs,   zfw,  zft 
    68       REAL(wp)                        :: zvmin, zamin, zamax  
    69       REAL(wp)                        :: zvtrp, zetrp 
    70       REAL(wp)                        :: zarea, zv_sill, zs_sill, zh_sill 
    71       REAL(wp), PARAMETER             :: zconv = 1.e-9 ! convert W to GW and kg to Mt 
     65      REAL(wp)        , INTENT(inout) ::   pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft 
     66      !! 
     67      REAL(wp) ::   zv, zs, zt, zfs, zfv, zft 
     68      REAL(wp) ::  zvmin, zamin, zamax  
     69      REAL(wp) ::  zvtrp, zetrp 
     70      REAL(wp) ::   zarea, zv_sill, zs_sill, zt_sill 
     71      REAL(wp), PARAMETER ::  zconv = 1.e-9 ! convert W to GW and kg to Mt 
    7272      !!---------------------------------------------------------------------- 
    7373      ! 
    7474      IF( icount == 0 ) THEN 
     75         !                          ! water flux 
     76         pdiag_fv = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:)     + wfx_sum(:,:)     + wfx_sni(:,:)     +                     & 
     77            &                    wfx_opw(:,:) + wfx_res(:,:)     + wfx_dyn(:,:)     + wfx_lam(:,:)     + wfx_ice_sub(:,:) +  & 
     78            &                    wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)    & 
     79            &                  ) * e1e2t(:,:) ) * zconv 
     80         ! 
    7581         !                          ! salt flux 
    76          zfs_b  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
    77             &                  sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)    & 
    78             &                ) *  e1e2t(:,:) ) * zconv  
    79          ! 
    80          !                          ! water flux 
    81          zfw_b  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:)     + wfx_sum(:,:)     + wfx_sni(:,:)     +                     & 
    82             &                  wfx_opw(:,:) + wfx_res(:,:)     + wfx_dyn(:,:)     + wfx_lam(:,:)     + wfx_ice_sub(:,:) +  & 
    83             &                  wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)    & 
    84             &                ) * e1e2t(:,:) ) * zconv 
     82         pdiag_fs = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
     83            &                    sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)    & 
     84            &                  ) *  e1e2t(:,:) ) * zconv  
    8585         ! 
    8686         !                          ! heat flux 
    87          zft_b = glob_sum(  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  &  
    88             &                - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   & 
    89             &                ) *  e1e2t(:,:) ) * zconv 
    90  
    91          zvi_b = glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t * zconv ) 
    92  
    93          zsmv_b = glob_sum( SUM( smv_i * rhoic            , dim=3 ) * e1e2t * zconv ) 
    94  
    95          zei_b = glob_sum( (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )     & 
    96             &                + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )   ) * e1e2t ) * zconv 
     87         pdiag_ft = glob_sum(  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  &  
     88            &                  - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   & 
     89            &                  ) *  e1e2t(:,:) ) * zconv 
     90 
     91         pdiag_v = glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t * zconv ) 
     92 
     93         pdiag_s = glob_sum( SUM( smv_i * rhoic            , dim=3 ) * e1e2t * zconv ) 
     94 
     95         pdiag_t = glob_sum( (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )     & 
     96            &                 + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )   ) * e1e2t ) * zconv 
    9797 
    9898      ELSEIF( icount == 1 ) THEN 
     99 
     100         ! water flux 
     101         zfv  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:)     + wfx_sum(:,:)     + wfx_sni(:,:)     +                     & 
     102            &                wfx_opw(:,:) + wfx_res(:,:)     + wfx_dyn(:,:)     + wfx_lam(:,:)     + wfx_ice_sub(:,:) +  & 
     103            &                wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)        & 
     104            &              ) * e1e2t(:,:) ) * zconv - pdiag_fv 
    99105 
    100106         ! salt flux 
    101107         zfs  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
    102108            &                sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)    &  
    103             &              ) * e1e2t(:,:) ) * zconv - zfs_b 
    104  
    105          ! water flux 
    106          zfw  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:)     + wfx_sum(:,:)     + wfx_sni(:,:)     +                     & 
    107             &                wfx_opw(:,:) + wfx_res(:,:)     + wfx_dyn(:,:)     + wfx_lam(:,:)     + wfx_ice_sub(:,:) +  & 
    108             &                wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)        & 
    109             &              ) * e1e2t(:,:) ) * zconv - zfw_b 
     109            &              ) * e1e2t(:,:) ) * zconv - pdiag_fs 
    110110 
    111111         ! heat flux 
    112112         zft  = glob_sum(  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  &  
    113113            &              - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   & 
    114             &              ) * e1e2t(:,:) ) * zconv - zft_b 
     114            &              ) * e1e2t(:,:) ) * zconv - pdiag_ft 
    115115  
    116116         ! outputs 
    117          zvi = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t  ) * zconv  & 
    118             &       - zvi_b ) * r1_rdtice - zfw ) * rday 
    119  
    120          zsmv = ( ( glob_sum( SUM( smv_i * rhoic            , dim=3 ) * e1e2t ) * zconv  & 
    121             &       - zsmv_b ) * r1_rdtice + zfs ) * rday 
    122  
    123          zei = ( glob_sum( (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )   & 
     117         zv = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t  ) * zconv  & 
     118            &     - pdiag_v ) * r1_rdtice - zfv ) * rday 
     119 
     120         zs = ( ( glob_sum( SUM( smv_i * rhoic            , dim=3 ) * e1e2t ) * zconv  & 
     121            &     - pdiag_s ) * r1_rdtice + zfs ) * rday 
     122 
     123         zt = ( glob_sum( (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )   & 
    124124            &                + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv   & 
    125             &   - zei_b ) * r1_rdtice + zft 
     125            &   - pdiag_t ) * r1_rdtice + zft 
    126126 
    127127         ! zvtrp and zetrp must be close to 0 if the advection scheme is conservative 
     
    137137         zv_sill = zarea * 2.5e-5 
    138138         zs_sill = zarea * 25.e-5 
    139          zh_sill = zarea * 10.e-5 
     139         zt_sill = zarea * 10.e-5 
    140140 
    141141         IF(lwp) THEN 
    142             IF ( ABS( zvi  ) > zv_sill ) WRITE(numout,*) 'violation volume [Mt/day]     (',cd_routine,') = ',zvi 
    143             IF ( ABS( zsmv ) > zs_sill ) WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zsmv 
    144             IF ( ABS( zei  ) > zh_sill ) WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',zei 
     142            IF ( ABS( zv   ) > zv_sill )   WRITE(numout,*) 'violation volume [Mt/day]     (',cd_routine,') = ',zv 
     143            IF ( ABS( zs   ) > zs_sill )   WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zs 
     144            IF ( ABS( zt   ) > zt_sill )   WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',zt 
    145145            IF ( ABS(zvtrp ) > zv_sill .AND. cd_routine == 'iceadv' ) THEN 
    146                                          WRITE(numout,*) 'violation vtrp [Mt/day]       (',cd_routine,') = ',zvtrp 
    147                                          WRITE(numout,*) 'violation etrp [GW]           (',cd_routine,') = ',zetrp 
    148             ENDIF 
    149             IF (     zvmin   < -epsi10 ) WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',zvmin 
    150             IF (     zamax   > MAX( rn_amax_n, rn_amax_s ) + epsi10 .AND. & 
    151                &                         cd_routine /= 'iceadv' .AND. cd_routine /= 'icerdgrft' ) THEN 
    152                                          WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax 
    153             IF (     zamax   > 1._wp   ) WRITE(numout,*) 'violation a_i>1               (',cd_routine,') = ',zamax 
    154             ENDIF 
    155             IF (      zamin  < -epsi10 ) WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin 
     146                                           WRITE(numout,*) 'violation vtrp [Mt/day]       (',cd_routine,') = ',zvtrp 
     147                                           WRITE(numout,*) 'violation etrp [GW]           (',cd_routine,') = ',zetrp 
     148            ENDIF 
     149            IF ( zvmin < -epsi10 )         WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',zvmin 
     150            IF ( zamax > MAX(rn_amax_n,rn_amax_s) + epsi10 .AND. cd_routine /= 'iceadv' .AND. cd_routine /= 'icerdgrft' )  & 
     151               &                           WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax 
     152            IF ( zamin < -epsi10 )         WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin 
    156153         ENDIF 
    157154         ! 
     
    169166      !! ** Method  : This is an online diagnostics which can be activated with ln_icediachk=true 
    170167      !!              It prints in ocean.output if there is a violation of conservation at each time-step 
    171       !!              The thresholds (zv_sill, zs_sill, zh_sill) which determine the violation are set to 
     168      !!              The thresholds (zv_sill, zs_sill, zt_sill) which determine the violation are set to 
    172169      !!              a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years. 
    173170      !!              For salt and heat thresholds, ice is considered to have a salinity of 10  
     
    176173      CHARACTER(len=*), INTENT(in)    :: cd_routine    ! name of the routine 
    177174      REAL(wp)                        :: zhfx, zsfx, zvfx 
    178       REAL(wp)                        :: zarea, zv_sill, zs_sill, zh_sill 
     175      REAL(wp)                        :: zarea, zv_sill, zs_sill, zt_sill 
    179176      REAL(wp), PARAMETER             :: zconv = 1.e-9 ! convert W to GW and kg to Mt 
    180177      !!---------------------------------------------------------------------- 
     178 
     179      ! water flux 
     180      zvfx  = glob_sum( ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) * zconv * rday 
     181 
     182      ! salt flux 
     183      zsfx  = glob_sum( ( sfx + diag_smvi ) * e1e2t ) * zconv * rday 
    181184 
    182185      ! heat flux 
     
    184187      !  &              - SUM( qevap_ice * a_i_b, dim=3 )                           & !!clem: I think this line must be commented (but need check) 
    185188         &              ) * e1e2t ) * zconv 
    186       ! salt flux 
    187       zsfx  = glob_sum( ( sfx + diag_smvi ) * e1e2t ) * zconv * rday 
    188       ! water flux 
    189       zvfx  = glob_sum( ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) * zconv * rday 
    190189 
    191190      ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice)  
     
    193192      zv_sill = zarea * 2.5e-5 
    194193      zs_sill = zarea * 25.e-5 
    195       zh_sill = zarea * 10.e-5 
     194      zt_sill = zarea * 10.e-5 
    196195 
    197196      IF(lwp) THEN 
    198          IF( ABS( zvfx ) > zv_sill )   WRITE(numout,*) 'violation vfx    [Mt/day]       (',cd_routine,')  = ',(zvfx) 
    199          IF( ABS( zsfx ) > zs_sill )   WRITE(numout,*) 'violation sfx    [psu*Mt/day]   (',cd_routine,')  = ',(zsfx) 
    200          IF( ABS( zhfx ) > zh_sill )   WRITE(numout,*) 'violation hfx    [GW]           (',cd_routine,')  = ',(zhfx) 
     197         IF( ABS( zvfx ) > zv_sill )   WRITE(numout,*) 'violation vfx  [Mt/day]       (',cd_routine,') = ',zvfx 
     198         IF( ABS( zsfx ) > zs_sill )   WRITE(numout,*) 'violation sfx  [psu*Mt/day]   (',cd_routine,') = ',zsfx 
     199         IF( ABS( zhfx ) > zt_sill )   WRITE(numout,*) 'violation hfx  [GW]           (',cd_routine,') = ',zhfx 
    201200      ENDIF 
    202201      ! 
     
    215214      INTEGER  ::   ialert_id         ! number of the current alert 
    216215      REAL(wp) ::   ztmelts           ! ice layer melting point 
    217       CHARACTER (len=30), DIMENSION(20)      ::   cl_alname   ! name of alert 
    218       INTEGER           , DIMENSION(20)      ::   inb_alp     ! number of alerts positive 
     216      CHARACTER (len=30), DIMENSION(20) ::   cl_alname   ! name of alert 
     217      INTEGER           , DIMENSION(20) ::   inb_alp     ! number of alerts positive 
    219218      !!------------------------------------------------------------------- 
    220219 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn.F90

    r8516 r8517  
    3434   PUBLIC   ice_dyn_init   ! called by icestp.F90 
    3535    
    36    INTEGER ::              nice_dyn   ! choice of the type of advection scheme 
     36   INTEGER ::              nice_dyn   ! choice of the type of dynamics 
    3737   !                                        ! associated indices: 
    38    INTEGER, PARAMETER ::   np_dynNO   = 0   ! no ice dynamics and ice advection 
    39    INTEGER, PARAMETER ::   np_dynFULL = 1   ! full ice dynamics  (rheology + advection + ridging/rafting + correction) 
    40    INTEGER, PARAMETER ::   np_dyn     = 2   ! no ridging/rafting (rheology + advection                   + correction) 
    41    INTEGER, PARAMETER ::   np_dynPURE = 3   ! pure dynamics      (rheology + advection)  
    42  
     38   INTEGER, PARAMETER ::   np_dynFULL    = 1   ! full ice dynamics               (rheology + advection + ridging/rafting + correction) 
     39   INTEGER, PARAMETER ::   np_dynRHGADV1 = 2   ! no ridging/rafting              (rheology + advection                   + correction) 
     40   INTEGER, PARAMETER ::   np_dynRHGADV2 = 3   ! pure dynamics                   (rheology + advection)  
     41   INTEGER, PARAMETER ::   np_dynADV     = 4   ! only advection w prescribed vel.(rn_uvice + advection) 
     42   ! 
     43   ! ** namelist (namdyn) ** 
     44   REAL(wp) ::   rn_uice          ! prescribed u-vel (case np_dynADV) 
     45   REAL(wp) ::   rn_vice          ! prescribed v-vel (case np_dynADV) 
     46    
    4347   !! * Substitutions 
    4448#  include "vectopt_loop_substitute.h90" 
     
    6165      INTEGER, INTENT(in) ::   kt     ! ice time step 
    6266      !! 
    63       INTEGER  ::   jl   ! dummy loop indices 
     67      INTEGER ::   ji, jj, jl         ! dummy loop indices 
     68      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhmax 
    6469      !!-------------------------------------------------------------------- 
    6570      ! 
     
    7277      ENDIF 
    7378 
    74       CALL ice_var_agg(1)           ! -- aggregate ice categories 
    7579      !                       
    76       IF( ln_landfast ) THEN        ! -- Landfast ice parameterization: define max bottom friction 
     80      IF( ln_landfast ) THEN            !-- Landfast ice parameterization: define max bottom friction 
    7781         tau_icebfr(:,:) = 0._wp 
    7882         DO jl = 1, jpl 
     
    8185      ENDIF 
    8286 
    83       SELECT CASE( nice_dyn )       ! -- Set which dynamics is running 
     87      zhmax(:,:,:) = ht_i_b(:,:,:)      !-- Record max of the surrounding 9-pts ice thick. (for CALL Hbig) 
     88      DO jl = 1, jpl 
     89         DO jj = 2, jpjm1 
     90            DO ji = 2, jpim1 
     91!!gm use of MAXVAL here is very probably less efficient than expending the 9 values 
     92               zhmax(ji,jj,jl) = MAX( epsi20, MAXVAL( ht_i_b(ji-1:ji+1,jj-1:jj+1,jl) ) ) 
     93            END DO 
     94         END DO 
     95      END DO 
     96      CALL lbc_lnk( zhmax(:,:,:), 'T', 1. ) 
     97      ! 
     98      ! 
     99      SELECT CASE( nice_dyn )           !-- Set which dynamics is running 
    84100 
    85101      CASE ( np_dynFULL )          !==  all dynamical processes  ==! 
    86          CALL ice_rhg   ( kt )          ! -- rheology   
    87          CALL ice_adv   ( kt )          ! -- advection of ice 
    88          CALL ice_rdgrft( kt )          ! -- ridging/rafting  
    89          CALL ice_cor   ( kt , 1 )      ! -- Corrections 
    90  
    91       CASE ( np_dyn )              !==  pure dynamics only ==!   (no ridging/rafting)   (nono cat. case 2) 
    92          CALL ice_rhg   ( kt )          ! -- rheology   
    93          CALL ice_adv   ( kt )          ! -- advection of ice 
    94          CALL ice_cor   ( kt , 1 )      ! -- Corrections 
    95  
    96       CASE ( np_dynPURE )          !==  pure dynamics only ==!   (nn_icedyn= 1 ) 
    97          CALL ice_rhg   ( kt )          ! -- rheology   
    98          CALL ice_adv   ( kt )          ! -- advection of ice 
    99  
    100       CASE ( np_dynNO )            !==  prescribed ice velocities ==!   (nn_icedyn= 0 ) 
     102         CALL ice_rhg   ( kt )                            ! -- rheology   
     103         CALL ice_adv   ( kt )   ;   CALL Hbig( zhmax )   ! -- advection of ice + correction on ice thickness 
     104         CALL ice_rdgrft( kt )                            ! -- ridging/rafting  
     105         CALL ice_cor   ( kt , 1 )                        ! -- Corrections 
     106 
     107      CASE ( np_dynRHGADV1 )       !==  no ridge/raft ==!   (mono cat. case 2) 
     108         CALL ice_rhg   ( kt )                            ! -- rheology   
     109         CALL ice_adv   ( kt )                            ! -- advection of ice 
     110         CALL Hpiling                                     ! -- simple pile-up (replaces ridging/rafting) 
     111         CALL ice_cor   ( kt , 1 )                        ! -- Corrections 
     112 
     113      CASE ( np_dynRHGADV2 )       !==  no ridge/raft & no corrections ==! 
     114         CALL ice_rhg   ( kt )                            ! -- rheology   
     115         CALL ice_adv   ( kt )                            ! -- advection of ice 
     116         CALL Hpiling                                     ! -- simple pile-up (replaces ridging/rafting) 
     117 
     118      CASE ( np_dynADV )           !==  pure advection ==!   (prescribed velocities) 
    101119         u_ice(:,:) = rn_uice * umask(:,:,1) 
    102120         v_ice(:,:) = rn_vice * vmask(:,:,1) 
    103121         !!CALL RANDOM_NUMBER(u_ice(:,:)) 
    104122         !!CALL RANDOM_NUMBER(v_ice(:,:)) 
     123         CALL ice_adv   ( kt )                            ! -- advection of ice 
    105124 
    106125      END SELECT 
     
    110129   END SUBROUTINE ice_dyn 
    111130 
     131   SUBROUTINE Hbig( phmax ) 
     132      !!------------------------------------------------------------------- 
     133      !!                  ***  ROUTINE Hbig  *** 
     134      !! 
     135      !! ** Purpose : Thickness correction in case advection scheme creates 
     136      !!              abnormally tick ice 
     137      !! 
     138      !! ** Method  : 1- check whether ice thickness resulting from advection is 
     139      !!                 larger than the surrounding 9-points before advection 
     140      !!                 and reduce it if a) divergence or b) convergence & at_i>0.8 
     141      !!              2- bound ice thickness with hi_max (99m) 
     142      !! 
     143      !! ** input   : Max thickness of the surrounding 9-points 
     144      !!------------------------------------------------------------------- 
     145      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phmax   ! max ice thick from surrounding 9-pts 
     146      ! 
     147      INTEGER  ::   ji, jj, jl         ! dummy loop indices 
     148      REAL(wp) ::   zh, zdv 
     149      !!------------------------------------------------------------------- 
     150      ! 
     151      CALL ice_var_zapsmall                       !-- zap small areas 
     152      ! 
     153      DO jl = 1, jpl 
     154         DO jj = 1, jpj 
     155            DO ji = 1, jpi 
     156               IF ( v_i(ji,jj,jl) > 0._wp ) THEN  !-- bound to hmax 
     157                  ! 
     158                  zh  = v_i (ji,jj,jl) / a_i(ji,jj,jl) 
     159                  zdv = v_i(ji,jj,jl) - v_i_b(ji,jj,jl)   
     160                  ! 
     161                  IF ( ( zdv >  0.0 .AND. zh > phmax(ji,jj,jl) .AND. at_i_b(ji,jj) < 0.80 ) .OR. & 
     162                     & ( zdv <= 0.0 .AND. zh > phmax(ji,jj,jl) ) ) THEN 
     163                     a_i (ji,jj,jl) = v_i(ji,jj,jl) / MIN( phmax(ji,jj,jl), hi_max(jpl) )   !-- bound ht_i to hi_max (99 m) 
     164                  ENDIF 
     165                  ! 
     166               ENDIF 
     167            END DO 
     168         END DO 
     169      END DO 
     170             
     171      IF ( nn_pnd_scheme > 0 ) THEN               !-- correct pond fraction to avoid a_ip > a_i 
     172         WHERE( a_ip(:,:,:) > a_i(:,:,:) )   a_ip(:,:,:) = a_i(:,:,:) 
     173      ENDIF 
     174      ! 
     175   END SUBROUTINE Hbig 
     176 
     177   SUBROUTINE Hpiling 
     178      !!------------------------------------------------------------------- 
     179      !!                  ***  ROUTINE Hpiling  *** 
     180      !! 
     181      !! ** Purpose : Simple conservative piling comparable with 1-cat models 
     182      !! 
     183      !! ** Method  : pile-up ice when no ridging/rafting 
     184      !! 
     185      !! ** input   : a_i 
     186      !!------------------------------------------------------------------- 
     187      INTEGER ::   jl         ! dummy loop indices 
     188      !!------------------------------------------------------------------- 
     189      ! 
     190      CALL ice_var_zapsmall                       !-- zap small areas 
     191      ! 
     192      at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 
     193      DO jl = 1, jpl 
     194         WHERE( at_i(:,:) > epsi20 ) 
     195            a_i(:,:,jl) = a_i(:,:,jl) * (  1._wp + MIN( rn_amax_2d(:,:) - at_i(:,:) , 0._wp ) / at_i(:,:)  ) 
     196         END WHERE 
     197      END DO 
     198      ! 
     199   END SUBROUTINE Hpiling 
    112200 
    113201   SUBROUTINE ice_dyn_init 
     
    123211      !! ** input   :   Namelist namice_dyn 
    124212      !!------------------------------------------------------------------- 
    125       INTEGER ::   ios   ! Local integer output status for namelist read 
    126       !! 
    127       NAMELIST/namice_dyn/ ln_icedyn  , nn_icedyn, rn_uice  , rn_vice ,    & 
     213      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read 
     214      !! 
     215      NAMELIST/namice_dyn/ ln_dynFULL, ln_dynRHGADV, ln_dynADV, rn_uice, rn_vice,  & 
    128216         &                 rn_ishlat  , rn_cio   ,                         & 
    129217         &                 ln_landfast, rn_gamma , rn_icebfr, rn_lfrelax 
     
    144232         WRITE(numout,*) '~~~~~~~~~~~~' 
    145233         WRITE(numout,*) '   Namelist namice_dyn' 
    146          WRITE(numout,*) '      Ice dynamics       (T) or not (F)                         ln_icedyn  = ', ln_icedyn 
    147          WRITE(numout,*) '         associated switch                                      nn_icedyn  = ', nn_icedyn 
    148          WRITE(numout,*) '            =2 all processes (default option)' 
    149          WRITE(numout,*) '            =1 advection only (no ridging/rafting)' 
    150          WRITE(numout,*) '            =0 advection only with prescribed velocity given by ' 
     234         WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr)  ln_dynFULL   = ', ln_dynFULL 
     235         WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                      ln_dynRHGADV = ', ln_dynRHGADV 
     236         WRITE(numout,*) '      Advection only         (rn_uvice + adv)                 ln_dynADV    = ', ln_dynADV 
     237         WRITE(numout,*) '           with prescribed velocity given by ' 
    151238         WRITE(numout,*) '               a uniform field               (u,v)_ice = (rn_uice,rn_vice) = (', rn_uice,',', rn_vice,')' 
    152239         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics        rn_ishlat     = ', rn_ishlat 
     
    158245      ENDIF 
    159246      !                             !== set the choice of ice dynamics ==! 
    160       SELECT CASE( nn_icedyn ) 
    161       CASE( 2 )                     
    162          IF( nn_monocat /= 2 ) THEN          !--- full dynamics (rheology + advection + ridging/rafting + correction) 
    163             nice_dyn = np_dynFULL 
    164          ELSE 
    165             nice_dyn = np_dyn                !--- dynamics without ridging/rafting 
    166          ENDIF 
    167       CASE( 1 )                              !--- dynamics without ridging/rafting and correction  
    168          nice_dyn = np_dynPURE 
    169       CASE( 0 )                              !--- prescribed ice velocities (from namelist)  
    170          nice_dyn = np_dynNO 
    171       END SELECT 
     247      ioptio = 0  
     248      !      !--- full dynamics                               (rheology + advection + ridging/rafting + correction) 
     249      IF( ln_dynFULL .AND. nn_monocat /= 2 ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynFULL      ;   ENDIF 
     250      !      !--- dynamics without ridging/rafting            (rheology + advection                   + correction) 
     251      IF( ln_dynFULL .AND. nn_monocat == 2 ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynRHGADV1   ;   ENDIF 
     252      !      !--- dynamics without ridging/rafting and corr   (rheology + advection) 
     253      IF( ln_dynRHGADV                     ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynRHGADV2   ;   ENDIF 
     254      !      !--- advection only with prescribed ice velocities (from namelist) 
     255      IF( ln_dynADV                        ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV       ;   ENDIF 
     256      ! 
     257      IF( ioptio /= 1 )    CALL ctl_stop( 'ice_dyn_init: one and only one ice dynamics option has to be defined ' ) 
    172258      ! 
    173259      !                                      !--- Lateral boundary conditions 
     
    177263      ELSEIF ( 2. < rn_ishlat                      ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  strong-slip' 
    178264      ENDIF 
    179       ! 
    180265      !                                      !--- NO Landfast ice : set to zero once for all 
    181266      IF( .NOT. ln_landfast )   tau_icebfr(:,:) = 0._wp  
    182267      ! 
    183       !                                      !--- simple conservative piling, comparable with LIM2 
    184       l_piling = nn_icedyn == 1 .OR. ( nn_monocat == 2  .AND.  jpl == 1 ) 
    185       ! 
    186268   END SUBROUTINE ice_dyn_init 
    187269 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90

    r8515 r8517  
    6767      REAL(wp) ::   zx3         
    6868      REAL(wp) ::   zslope          ! used to compute local thermodynamic "speeds" 
    69       REAL(wp) ::   zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b ! conservation check 
    7069 
    7170      INTEGER , DIMENSION(jpij)       ::   idxice2         ! compute remapping or not 
     
    8180      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_rem: remapping ice thickness distribution'  
    8281 
    83       IF( ln_icediachk ) CALL ice_cons_hsm(0, 'iceitd_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     82      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
    8483 
    8584      !----------------------------------------------------------------------------------------------- 
     
    311310      ENDIF 
    312311      ! 
    313       IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     312      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
    314313      ! 
    315314   END SUBROUTINE ice_itd_rem 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90

    r8515 r8517  
    127127      ! 
    128128      INTEGER, PARAMETER ::   nitermax = 20     
    129       ! 
    130       REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    131129      !!----------------------------------------------------------------------------- 
    132       IF( nn_timing == 1 )  CALL timing_start('icerdgrft') 
     130      ! controls 
     131      IF( nn_timing == 1 )   CALL timing_start('icerdgrft')                                                             ! timing 
     132      IF( ln_icediachk   )   CALL ice_cons_hsm(0, 'icerdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    133133 
    134134      IF( kt == nit000 ) THEN 
     
    140140         ! 
    141141      ENDIF       
    142       !                    ! conservation test 
    143       IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icerdgrft', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    144142 
    145143      !-----------------------------------------------------------------------------! 
     
    284282      CALL ice_var_agg( 1 )  
    285283 
    286       !-----------------------------------------------------------------------------! 
    287       ! control prints 
    288       !-----------------------------------------------------------------------------! 
    289       !                    ! conservation test 
    290       IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icerdgrft', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    291  
    292       !                    ! control prints 
    293       IF( ln_ctl       )   CALL ice_prt3D( 'icerdgrft' ) 
    294       ! 
    295       IF( nn_timing == 1 )  CALL timing_stop('icerdgrft') 
     284      ! controls 
     285      IF( ln_icediachk   )   CALL ice_cons_hsm(1, 'icerdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
     286      IF( ln_ctl         )   CALL ice_prt3D   ('icerdgrft')                                                             ! prints 
     287      IF( nn_timing == 1 )   CALL timing_stop ('icerdgrft')                                                             ! timing 
    296288      ! 
    297289   END SUBROUTINE ice_rdgrft 
     
    795787      ELSEIF( ln_str_H79 ) THEN      ! Ice strength => Hibler (1979) method             ! 
    796788      !                              !--------------------------------------------------! 
    797          strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) )  ) * tmask(:,:,1) 
     789         strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) ) 
    798790         ! 
    799791         ismooth = 1 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg.F90

    r8516 r8517  
    3232   PUBLIC   ice_rhg        ! called by icestp.F90 
    3333   PUBLIC   ice_rhg_init   ! called by icestp.F90 
     34 
     35   INTEGER ::              nice_rhg   ! choice of the type of rheology 
     36   !                                        ! associated indices: 
     37   INTEGER, PARAMETER ::   np_rhgEVP = 1   ! EVP rheology 
     38!! INTEGER, PARAMETER ::   np_rhgEAP = 2   ! EAP rheology 
    3439    
    3540   !! * Substitutions 
     
    5358      !!-------------------------------------------------------------------- 
    5459      INTEGER, INTENT(in) ::   kt     ! ice time step 
    55       !! 
    56       REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    5760      !!-------------------------------------------------------------------- 
    58       ! 
    59       IF( nn_timing == 1 )  CALL timing_start('icerhg') 
     61      ! controls 
     62      IF( nn_timing == 1 )   CALL timing_start('icerhg')                                                             ! timing 
     63      IF( ln_icediachk   )   CALL ice_cons_hsm(0, 'icerhg', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    6064      ! 
    6165      IF( kt == nit000 .AND. lwp ) THEN 
     
    6468         WRITE(numout,*)'~~~~~~~' 
    6569      ENDIF 
    66       !                             ! -- conservation test 
    67       IF( ln_icediachk   )   CALL ice_cons_hsm(0, 'icerhg', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    6870 
    69       ! ----------------------- 
    70       ! Rheology (ice dynamics) 
    71       ! -----------------------    
    72       CALL ice_rhg_evp( kt, stress1_i, stress2_i, stress12_i, u_ice, v_ice, shear_i, divu_i, delta_i ) 
     71      ! -------- 
     72      ! Rheology 
     73      ! --------    
     74      SELECT CASE( nice_rhg ) 
     75      !                                !------------------------! 
     76      CASE( np_rhgEVP )                ! Elasto-Viscous-Plastic ! 
     77         !                             !------------------------! 
     78         CALL ice_rhg_evp( kt, stress1_i, stress2_i, stress12_i, u_ice, v_ice, shear_i, divu_i, delta_i ) 
     79 
     80      END SELECT 
    7381      ! 
    74       !                             !- conservation test 
    75       IF( ln_icediachk   )   CALL ice_cons_hsm(1, 'icerhg', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    76       IF( ln_ctl         )   CALL ice_prt3D  ('icerhg')   !- Control prints 
    77       IF( nn_timing == 1 )   CALL timing_stop('icerhg')   !- timing 
     82      ! controls 
     83      IF( ln_icediachk   )   CALL ice_cons_hsm(1, 'icerhg', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
     84      IF( ln_ctl         )   CALL ice_prt3D   ('icerhg')                                                             ! prints 
     85      IF( nn_timing == 1 )   CALL timing_stop ('icerhg')                                                             ! timing 
    7886      ! 
    7987   END SUBROUTINE ice_rhg 
     
    92100      !! ** input   :   Namelist namice_rhg 
    93101      !!------------------------------------------------------------------- 
    94       INTEGER ::   ios   ! Local integer output status for namelist read 
     102      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read 
    95103      !! 
    96104      NAMELIST/namice_rhg/  ln_rhg_EVP, rn_creepl, rn_ecc , nn_nevp, rn_relast 
     
    118126      ENDIF 
    119127      ! 
     128      !                             !== set the choice of ice advection ==! 
     129      ioptio = 0  
     130      IF( ln_rhg_EVP ) THEN   ;   ioptio = ioptio + 1   ;   nice_rhg = np_rhgEVP    ;   ENDIF 
     131!!    IF( ln_rhg_EAP ) THEN   ;   ioptio = ioptio + 1   ;   nice_rhg = np_rhgEAP    ;   ENDIF 
     132      IF( ioptio /= 1 )   CALL ctl_stop( 'ice_rhg_init: choose one and only one ice rheology' ) 
     133 
     134     !                              ! allocate tke arrays 
     135!!clem example      IF( zdf_tke_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tke_init : unable to allocate arrays' ) 
     136      ! 
    120137   END SUBROUTINE ice_rhg_init 
    121138 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90

    r8515 r8517  
    2323   USE sbc_oce , ONLY : ln_ice_embd, nn_fsbc, ssh_m 
    2424   USE sbc_ice , ONLY : utau_ice, vtau_ice, snwice_mass, snwice_mass_b 
    25    USE ice            ! ice variables 
    26    USE icerdgrft      ! ice strength 
     25   USE ice            ! sea-ice: ice variables 
     26   USE icerdgrft      ! sea-ice: ice strength 
    2727   ! 
    2828   USE lbclnk         ! Lateral Boundary Condition / MPP link 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90

    r8516 r8517  
    258258      CALL ice_itd_init                ! ice thickness distribution initialization 
    259259      ! 
    260       CALL ice_dyn_init                ! set ice dynamics parameters 
    261       ! 
    262       CALL ice_rdgrft_init             ! set ice ridging/rafting parameters 
    263       ! 
    264       CALL ice_rhg_init                ! set ice rheology parameters 
    265       ! 
    266       CALL ice_adv_init                ! set ice advection parameters 
    267       ! 
    268       CALL ice_thd_init                ! set ice thermodynics parameters 
    269       ! 
    270       CALL ice_thd_sal_init            ! set ice salinity parameters 
    271         
     260      IF( ln_icedyn )   THEN 
     261         CALL ice_dyn_init             ! set ice dynamics parameters 
     262         CALL ice_rdgrft_init          ! set ice ridging/rafting parameters 
     263         CALL ice_rhg_init             ! set ice rheology parameters 
     264         CALL ice_adv_init             ! set ice advection parameters 
     265      ENDIF 
     266 
     267      IF( ln_icethd ) THEN 
     268         CALL ice_thd_init             ! set ice thermodynics parameters 
     269         CALL ice_thd_sal_init         ! set ice salinity parameters 
     270      ENDIF 
     271    
    272272      ! MV MP 2016 
    273273      CALL lim_mp_init                 ! set melt ponds parameters 
     
    283283      CALL ice_var_glo2eqv 
    284284      ! 
    285       CALL ice_update_init                 ! ice surface boundary condition 
    286       ! 
    287       CALL ice_alb_init                    ! ice surface albedo 
    288       ! 
    289       CALL ice_dia_init                    ! initialization for diags 
    290       ! 
    291       fr_i  (:,:)   = at_i(:,:)         ! initialisation of sea-ice fraction 
    292       tn_ice(:,:,:) = t_su(:,:,:)       ! initialisation of surface temp for coupled simu 
    293       ! 
     285      CALL ice_update_init             ! ice surface boundary condition 
     286      ! 
     287      CALL ice_alb_init                ! ice surface albedo 
     288      ! 
     289      CALL ice_dia_init                ! initialization for diags 
     290      ! 
     291      fr_i  (:,:)   = at_i(:,:)        ! initialisation of sea-ice fraction 
     292      tn_ice(:,:,:) = t_su(:,:,:)      ! initialisation of surface temp for coupled simu 
     293      ! 
     294      !                                ! set max concentration in both hemispheres 
    294295      WHERE( gphit(:,:) > 0._wp )   ;   rn_amax_2d(:,:) = rn_amax_n  ! NH 
    295296      ELSEWHERE                     ;   rn_amax_2d(:,:) = rn_amax_s  ! SH 
     
    306307      !! 
    307308      !! ** Method  :   Read the namice_run namelist and check the parameter 
    308       !!              values called at the first timestep (nit000) 
     309      !!                values called at the first timestep (nit000) 
    309310      !! 
    310311      !! ** input   :   Namelist namice_run 
    311312      !!------------------------------------------------------------------- 
    312313      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    313       NAMELIST/namice_run/ jpl, nlay_i, nlay_s, nn_monocat, rn_amax_n, rn_amax_s, cn_icerst_in, cn_icerst_indir,   & 
    314          &                cn_icerst_out, cn_icerst_outdir 
     314      NAMELIST/namice_run/ jpl, nlay_i, nlay_s, nn_monocat, ln_icedyn, ln_icethd, rn_amax_n, rn_amax_s,  & 
     315         &                 cn_icerst_in, cn_icerst_indir, cn_icerst_out, cn_icerst_outdir 
    315316      !!------------------------------------------------------------------- 
    316317      ! 
     
    333334         WRITE(numout,*) '      number of snow layers                                  nlay_s = ', nlay_s 
    334335         WRITE(numout,*) '      virtual ITD mono-category param (1-4) or not (0)   nn_monocat = ', nn_monocat 
     336         WRITE(numout,*) '      Ice dynamics       (T) or not (F)                  ln_icedyn  = ', ln_icedyn 
     337         WRITE(numout,*) '      Ice thermodynamics (T) or not (F)                  ln_icethd  = ', ln_icethd 
    335338         WRITE(numout,*) '      maximum ice concentration for NH                              = ', rn_amax_n  
    336339         WRITE(numout,*) '      maximum ice concentration for SH                              = ', rn_amax_s 
    337340      ENDIF 
    338341      ! 
    339       !                                         !--- check consistency 
     342      !                                        !--- check consistency 
    340343      IF ( jpl > 1 .AND. nn_monocat == 1 ) THEN 
    341344         nn_monocat = 0 
     
    349352      IF( ln_bdy .AND. ln_icediachk )   CALL ctl_warn('ice_run_init: online conservation check does not work with BDY') 
    350353      ! 
    351       rdt_ice   = REAL(nn_fsbc) * rdt           !--- sea-ice timestep and inverse 
     354      rdt_ice   = REAL(nn_fsbc) * rdt          !--- sea-ice timestep and inverse 
    352355      r1_rdtice = 1._wp / rdt_ice 
    353       IF( lwp ) WRITE(numout,*) '   ice timestep rdt_ice  = ', rdt_ice 
    354       ! 
    355       r1_nlay_i = 1._wp / REAL( nlay_i, wp )    !--- inverse of nlay_i and nlay_s 
     356      IF( lwp ) WRITE(numout,*) '   ice timestep rdt_ice = ', rdt_ice 
     357      ! 
     358      r1_nlay_i = 1._wp / REAL( nlay_i, wp )   !--- inverse of nlay_i and nlay_s 
    356359      r1_nlay_s = 1._wp / REAL( nlay_s, wp ) 
    357360      ! 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90

    r8515 r8517  
    8181      INTEGER  :: ji, jj, jk, jl   ! dummy loop indices 
    8282      REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg 
    83       REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b ! conservation check 
    8483      REAL(wp), PARAMETER :: zfric_umin = 0._wp           ! lower bound for the friction velocity (cice value=5.e-04) 
    8584      REAL(wp), PARAMETER :: zch        = 0.0057_wp       ! heat transfer coefficient 
     
    8786      ! 
    8887      !!------------------------------------------------------------------- 
    89  
    90       IF( nn_timing == 1 )   CALL timing_start('icethd') 
     88      ! controls 
     89      IF( nn_timing == 1 )   CALL timing_start('icethd')                                                             ! timing 
     90      IF( ln_icediachk   )   CALL ice_cons_hsm(0, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    9191 
    9292      IF( kt == nit000 .AND. lwp ) THEN 
     
    9696      ENDIF 
    9797       
    98       ! conservation test 
    99       IF( ln_icediachk )   CALL ice_cons_hsm( 0, 'icethd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b ) 
    100  
    10198      CALL ice_var_glo2eqv 
    10299 
     
    257254      oa_i(:,:,:) = o_i(:,:,:) * a_i(:,:,:) 
    258255 
    259       IF( ln_icediachk )   CALL ice_cons_hsm( 1, 'icethd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b ) 
     256      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
    260257      ! 
    261258                           CALL ice_var_zapsmall           ! --- remove very small ice concentration (<1e-10) --- ! 
     
    266263      IF( ln_icedO )       CALL ice_thd_lac                ! --- frazil ice growing in leads --- ! 
    267264      ! 
    268       IF( ln_icectl )      CALL ice_prt( kt, iiceprt, jiceprt, 1, ' - ice thermodyn. - ' )   ! control print 
    269       IF( ln_ctl )         CALL ice_prt3D( 'icethd' )      ! Control print 
    270       ! 
    271       IF( nn_timing == 1 )  CALL timing_stop('icethd') 
     265      ! controls 
     266      IF( ln_icectl      )   CALL ice_prt    (kt, iiceprt, jiceprt, 1, ' - ice thermodyn. - ') ! prints 
     267      IF( ln_ctl         )   CALL ice_prt3D  ('icethd')                                        ! prints 
     268      IF( nn_timing == 1 )   CALL timing_stop('icethd')                                        ! timing 
    272269      ! 
    273270   END SUBROUTINE ice_thd  
     
    539536      INTEGER  ::   ios   ! Local integer output status for namelist read 
    540537      !! 
    541       NAMELIST/namice_thd/ ln_icethd, rn_kappa_i, ln_cndi_U64, ln_cndi_P07, ln_dqns_i, rn_cnd_s,   & 
     538      NAMELIST/namice_thd/ rn_kappa_i, ln_cndi_U64, ln_cndi_P07, ln_dqns_i, rn_cnd_s,   & 
    542539         &                 ln_icedH, rn_blow_s,                                                    & 
    543540         &                 ln_icedA, rn_beta, rn_dmin,                                             & 
     
    560557         WRITE(numout,*) '~~~~~~~~~~~~~' 
    561558         WRITE(numout,*) '   Namelist namice_thd' 
    562          WRITE(numout,*) '      Ice thermodynamics (T) or not (F)                       ln_icethd    = ', ln_icethd 
    563559         WRITE(numout,*) '   -- icethd_dif --' 
    564560         WRITE(numout,*) '      extinction radiation parameter in sea ice               rn_kappa_i   = ', rn_kappa_i 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_lac.F90

    r8514 r8517  
    8383      
    8484      REAL(wp) ::   zv_newfra 
    85       REAL(wp) ::   zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b ! conservation check 
    8685   
    8786      INTEGER , DIMENSION(jpij) ::   jcat        ! indexes of categories where new ice grows 
     
    113112      !!-----------------------------------------------------------------------! 
    114113 
    115       IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icethd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     114      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icethd_lac', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
    116115 
    117116      CALL ice_var_agg(1) 
     
    479478      ENDIF ! nidx > 0 
    480479      ! 
    481       IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icethd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     480      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd_lac', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
    482481      ! 
    483482   END SUBROUTINE ice_thd_lac 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90

    r8514 r8517  
    104104      ! END MP 2016 
    105105 
    106       DO jj = 1, jpj                         ! open water fraction 
    107          DO ji = 1, jpi 
    108             ato_i(ji,jj) = MAX( 1._wp - at_i(ji,jj), 0._wp )    
    109          END DO 
    110       END DO 
    111 !!gm I think this should do the work : 
    112 !      ato_i(:,:) = MAX( 1._wp - at_i(:,:), 0._wp )   
    113 !!gm end 
     106      ato_i(:,:) = 1._wp - at_i(:,:)    ! open water fraction   
    114107 
    115108      IF( kn > 1 ) THEN 
     
    543536 
    544537      ! to be sure that at_i is the sum of a_i(jl) 
    545       at_i (:,:) = a_i(:,:,1) 
    546       vt_i (:,:) = v_i(:,:,1) 
    547       DO jl = 2, jpl 
    548          at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 
    549          vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 
    550       END DO 
    551  
    552       ! open water = 1 if at_i=0 (no re-calculation of ato_i here) 
    553       DO jj = 1, jpj 
    554          DO ji = 1, jpi 
    555             rswitch      = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) ) 
    556             ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj) 
    557          END DO 
    558       END DO 
     538      at_i (:,:) = SUM( a_i(:,:,:), dim=3 ) 
     539      vt_i (:,:) = SUM( v_i(:,:,:), dim=3 ) 
     540 
     541      ! open water = 1 if at_i=0 
     542      WHERE( at_i(:,:) == 0._wp )   ato_i(:,:) = 1._wp 
    559543      ! 
    560544   END SUBROUTINE ice_var_zapsmall 
Note: See TracChangeset for help on using the changeset viewer.