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 8486 – NEMO

Changeset 8486


Ignore:
Timestamp:
2017-09-01T15:49:35+02:00 (7 years ago)
Author:
clem
Message:

changes in style - part1 - (now the code looks better txs to Gurvan's comments)

Location:
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3
Files:
28 edited

Legend:

Unmodified
Added
Removed

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r8426 r8486  
    1717   PRIVATE 
    1818 
    19    PUBLIC    ice_alloc  !  Called in ice_init 
     19   PUBLIC    ice_alloc  ! called by icestp.F90 
    2020 
    2121   !!====================================================================== 
    … …  
    288288   REAL(wp), PUBLIC ::   r1_nlay_s        !: 1 / nlay_s  
    289289   REAL(wp), PUBLIC ::   rswitch          !: switch for the presence of ice (1) or not (0) 
    290    REAL(wp), PUBLIC, PARAMETER ::   epsi06   = 1.e-06_wp  !: small number  
    291    REAL(wp), PUBLIC, PARAMETER ::   epsi10   = 1.e-10_wp  !: small number  
    292    REAL(wp), PUBLIC, PARAMETER ::   epsi20   = 1.e-20_wp  !: small number  
     290   REAL(wp), PUBLIC, PARAMETER ::   epsi06 = 1.e-06_wp  !: small number  
     291   REAL(wp), PUBLIC, PARAMETER ::   epsi10 = 1.e-10_wp  !: small number  
     292   REAL(wp), PUBLIC, PARAMETER ::   epsi20 = 1.e-20_wp  !: small number  
     293   ! 
     294   LOGICAL , PUBLIC ::   l_piling         !: =T simple conservative piling, comparable with LIM2 
    293295 
    294296   !                                     !!** define arrays 
    295    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   u_oce, v_oce !: surface ocean velocity used in ice dynamics 
     297   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   u_oce,v_oce !: surface ocean velocity used in ice dynamics 
    296298   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hicol       !: ice collection thickness accreted in leads 
    297299   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   strength    !: ice strength 
    … …  
    499501   ! 
    500502   !!---------------------------------------------------------------------- 
    501    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 
     503   !! NEMO/ICE 4.0 , NEMO Consortium (2017) 
    502504   !! $Id$ 
    503505   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1D.F90

    r8420 r8486  
    88#if defined key_lim3 
    99   !!---------------------------------------------------------------------- 
    10    !!   'key_lim3'                                      LIM3 sea-ice model 
    11    !!---------------------------------------------------------------------- 
     10   !!   'key_lim3'                                       LIM3 sea-ice model 
     11   !!---------------------------------------------------------------------- 
     12   USE ice     , ONLY :   nlay_i, nlay_s, jpl   ! number of ice/snow layers and categories 
     13   ! 
    1214   USE in_out_manager ! I/O manager 
    1315   USE lib_mpp        ! MPP library 
    14    USE ice, ONLY :   nlay_i, nlay_s, jpl 
    1516 
    1617   IMPLICIT NONE 
    1718   PRIVATE 
    1819 
    19    PUBLIC ice1D_alloc ! Routine called by nemogcm.F90 
     20   PUBLIC   ice1D_alloc   ! called by icestp.F90 
    2021 
    2122   !!---------------------- 
    … …  
    2627   !: are the variables corresponding to 2d vectors 
    2728 
     29!! please, DOCTOR naming convention   starting by i means LOCAL integer 
     30!!         ===>>>  rename idxice  as nidice or nidthd, or nidx_thd or midx ... ? 
    2831   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   idxice !: selected points for ice thermo 
    2932   INTEGER , PUBLIC                                  ::   nidx   !  number of selected points 
    … …  
    151154    
    152155   !!---------------------------------------------------------------------- 
    153    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     156   !! NEMO/ICE 4.0 , NEMO Consortium (2017) 
    154157   !! $Id: ice1D.F90 8379 2017-07-26 15:35:49Z clem $ 
    155158   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90

    r8426 r8486  
    1010#if defined key_lim3 
    1111   !!---------------------------------------------------------------------- 
    12    !!   'key_lim3'                                      LIM3 sea-ice model 
    13    !!---------------------------------------------------------------------- 
    14    !!   ice_adv      : advection/diffusion process of sea ice 
     12   !!   'key_lim3'                                       LIM3 sea-ice model 
     13   !!---------------------------------------------------------------------- 
     14   !!   ice_adv       : advection/diffusion process of sea ice 
    1515   !!---------------------------------------------------------------------- 
    1616   USE phycst         ! physical constant 
    1717   USE dom_oce        ! ocean domain 
    18    USE sbc_oce , ONLY : nn_fsbc 
    19    USE ice            ! ice variables 
    20    USE icevar         !  
    21    USE iceadv_prather ! advection scheme (Prather) 
    22    USE iceadv_umx     ! advection scheme (ultimate-macho) 
     18   USE sbc_oce , ONLY : nn_fsbc   ! frequency of sea-ice call 
     19   USE ice            ! sea-ice: variables 
     20   USE icevar         ! sea-ice: operations 
     21   USE iceadv_prather ! sea-ice: advection scheme (Prather) 
     22   USE iceadv_umx     ! sea-ice: advection scheme (ultimate-macho) 
     23   USE icectl         ! sea-ice: control prints 
    2324   ! 
    2425   USE in_out_manager ! I/O manager 
    2526   USE lbclnk         ! lateral boundary conditions -- MPP exchanges 
    2627   USE lib_mpp        ! MPP library 
    27    USE wrk_nemo       ! work arrays 
    2828   USE prtctl         ! Print control 
    2929   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3030   USE timing         ! Timing 
    31    USE icectl         ! control prints 
    3231 
    3332   IMPLICIT NONE 
    … …  
    3635   PUBLIC   ice_adv    ! called by icestp 
    3736 
    38    INTEGER  ::   ncfl                 ! number of ice time step with CFL>1/2   
     37   INTEGER ::   ncfl   ! number of ice time step with CFL>1/2   
    3938 
    4039   !! * Substitution 
    4140#  include "vectopt_loop_substitute.h90" 
    4241   !!---------------------------------------------------------------------- 
    43    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     42   !! NEMO/ICE 4.0 , NEMO Consortium (2017) 
    4443   !! $Id: iceadv.F90 8373 2017-07-25 17:44:54Z clem $ 
    4544   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    4847 
    4948   SUBROUTINE ice_adv( kt )  
    50       !!------------------------------------------------------------------- 
     49      !!---------------------------------------------------------------------- 
    5150      !!                   ***  ROUTINE ice_adv *** 
    5251      !!                     
    … …  
    6059      !! 
    6160      !! ** action : 
    62       !!--------------------------------------------------------------------- 
     61      !!---------------------------------------------------------------------- 
    6362      INTEGER, INTENT(in) ::   kt   ! number of iteration 
    6463      ! 
    … …  
    7372      REAL(wp), DIMENSION(jpi,jpj,jpl)       ::   zhimax, zviold, zvsold 
    7473      ! --- ultimate macho only --- ! 
    75       REAL(wp)                               ::   zdt 
    76       REAL(wp), POINTER, DIMENSION(:,:)      ::   zudy, zvdx, zcu_box, zcv_box 
     74      REAL(wp) ::   zdt 
     75      REAL(wp), ALLOCATABLE, DIMENSION(:,:)      ::   zudy, zvdx, zcu_box, zcv_box 
    7776      ! --- prather only --- ! 
    78       REAL(wp), POINTER, DIMENSION(:,:)      ::   zarea 
    79       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0opw 
    80       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0ice, z0snw, z0ai, z0es , z0smi , z0oi 
     77      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   zarea 
     78      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0opw 
     79      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0ice, z0snw, z0ai, z0es , z0smi , z0oi 
    8180      ! MV MP 2016 
    82       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0ap , z0vp 
     81      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0ap , z0vp 
    8382      REAL(wp) ::   za_old 
    8483      ! END MV MP 2016 
    85       REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::   z0ei 
    86       !! 
     84      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   z0ei 
    8785      !!--------------------------------------------------------------------- 
    8886      IF( nn_timing == 1 )  CALL timing_start('iceadv') 
    … …  
    9088      IF( kt == nit000 .AND. lwp ) THEN 
    9189         WRITE(numout,*)'' 
    92          WRITE(numout,*)'iceadv' 
    93          WRITE(numout,*)'~~~~~~' 
     90         WRITE(numout,*)'iceadv : sea-ice advection' 
     91         WRITE(numout,*)'~~~~~~~' 
    9492         ncfl = 0                ! nb of time step with CFL > 1/2 
    9593      ENDIF 
    … …  
    127125         DO jj = 2, jpjm1 
    128126            DO ji = 2, jpim1 
     127!!gm use of MAXVAL here is very probably less efficient than expending the 9 values 
    129128               zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) + ht_s(ji-1:ji+1,jj-1:jj+1,jl) ) 
    130129            END DO 
    … …  
    156155                       !=============================! 
    157156       
    158          CALL wrk_alloc( jpi,jpj, zudy, zvdx, zcu_box, zcv_box ) 
     157         ALLOCATE( zudy(jpi,jpj) , zvdx(jpi,jpj) , zcu_box(jpi,jpj) , zcv_box(jpi,jpj) ) 
    159158       
    160159         IF( kt == nit000 .AND. lwp ) THEN 
    … …  
    212211         END DO 
    213212         ! 
    214          CALL wrk_dealloc( jpi,jpj, zudy, zvdx, zcu_box, zcv_box ) 
     213         DEALLOCATE( zudy, zvdx, zcu_box, zcv_box ) 
    215214          
    216215                       !=============================! 
    … …  
    218217                       !=============================! 
    219218 
    220          CALL wrk_alloc( jpi,jpj,            zarea ) 
    221          CALL wrk_alloc( jpi,jpj,1,          z0opw ) 
    222          CALL wrk_alloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi ) 
    223          CALL wrk_alloc( jpi,jpj,jpl,        z0ap , z0vp ) 
    224          CALL wrk_alloc( jpi,jpj,nlay_i,jpl, z0ei ) 
     219         ALLOCATE( zarea(jpi,jpj)     , z0opw(jpi,jpj, 1 ) ,                                           & 
     220            &      z0ice(jpi,jpj,jpl) , z0snw(jpi,jpj,jpl) , z0ai(jpi,jpj,jpl) , z0es(jpi,jpj,jpl) ,   & 
     221            &      z0smi(jpi,jpj,jpl) , z0oi (jpi,jpj,jpl) , z0ap(jpi,jpj,jpl) , z0vp(jpi,jpj,jpl) ,   & 
     222            &      z0ei (jpi,jpj,nlay_i,jpl)                                                           ) 
    225223          
    226224         IF( kt == nit000 .AND. lwp ) THEN 
    … …  
    237235         z0opw(:,:,1) = ato_i(:,:) * e1e2t(:,:)             ! Open water area  
    238236         DO jl = 1, jpl 
    239             z0snw (:,:,jl) = v_s  (:,:,  jl) * e1e2t(:,:)  ! Snow volume 
    240             z0ice(:,:,jl)   = v_i  (:,:,  jl) * e1e2t(:,:)  ! Ice  volume 
    241             z0ai  (:,:,jl) = a_i  (:,:,  jl) * e1e2t(:,:)  ! Ice area 
    242             z0smi (:,:,jl) = smv_i(:,:,  jl) * e1e2t(:,:)  ! Salt content 
    243             z0oi (:,:,jl)   = oa_i (:,:,  jl) * e1e2t(:,:)  ! Age content 
    244             z0es (:,:,jl)   = e_s  (:,:,1,jl) * e1e2t(:,:)  ! Snow heat content 
     237            z0snw(:,:,jl) = v_s  (:,:,  jl) * e1e2t(:,:)  ! Snow volume 
     238            z0ice(:,:,jl) = v_i  (:,:,  jl) * e1e2t(:,:)  ! Ice  volume 
     239            z0ai (:,:,jl) = a_i  (:,:,  jl) * e1e2t(:,:)  ! Ice area 
     240            z0smi(:,:,jl) = smv_i(:,:,  jl) * e1e2t(:,:)  ! Salt content 
     241            z0oi (:,:,jl) = oa_i (:,:,  jl) * e1e2t(:,:)  ! Age content 
     242            z0es (:,:,jl) = e_s  (:,:,1,jl) * e1e2t(:,:)  ! Snow heat content 
    245243            DO jk = 1, nlay_i 
    246                z0ei  (:,:,jk,jl) = e_i  (:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content 
     244               z0ei(:,:,jk,jl) = e_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content 
    247245            END DO 
    248246            ! MV MP 2016 
    249247            IF ( nn_pnd_scheme > 0 ) THEN 
    250                z0ap  (:,:,jl)  = a_ip (:,:,jl) * e1e2t(:,:) ! Melt pond fraction 
    251                z0vp  (:,:,jl)  = v_ip (:,:,jl) * e1e2t(:,:) ! Melt pond volume 
     248               z0ap(:,:,jl)  = a_ip(:,:,jl) * e1e2t(:,:) ! Melt pond fraction 
     249               z0vp(:,:,jl)  = v_ip(:,:,jl) * e1e2t(:,:) ! Melt pond volume 
    252250            ENDIF 
    253251            ! END MV MP 2016 
    254  
    255          END DO 
    256  
     252         END DO 
    257253 
    258254         IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==! 
    … …  
    383379            ! END MV MP 2016 
    384380         END DO 
    385  
     381         ! 
    386382         at_i(:,:) = a_i(:,:,1)      ! total ice fraction 
    387383         DO jl = 2, jpl 
    388384            at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 
    389385         END DO 
    390           
    391          CALL wrk_dealloc( jpi,jpj,            zarea ) 
    392          CALL wrk_dealloc( jpi,jpj,1,          z0opw ) 
    393          CALL wrk_dealloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi ) 
    394          ! MV MP 2016 
    395          CALL wrk_dealloc( jpi,jpj,jpl,        z0ap , z0vp ) 
    396          ! END MV MP 2016 
    397          CALL wrk_dealloc( jpi,jpj,nlay_i,jpl, z0ei ) 
    398  
     386         ! 
     387         DEALLOCATE( zarea , z0opw , z0ice, z0snw , z0ai , z0es , z0smi , z0oi , z0ap , z0vp , z0ei ) 
     388         ! 
    399389      END SELECT 
    400390       
    … …  
    411401       
    412402      IF( nn_limdyn == 2) THEN 
    413  
    414          ! zap small areas 
    415          CALL ice_var_zapsmall 
    416             
    417          !--- Thickness correction in case too high --- ! 
    418          DO jl = 1, jpl 
     403         ! 
     404         CALL ice_var_zapsmall      !--- zap small areas ---! 
     405         ! 
     406         DO jl = 1, jpl             !--- Thickness correction in case too high --- ! 
    419407            DO jj = 1, jpj 
    420408               DO ji = 1, jpi 
    421                    
     409                  !  
    422410                  IF ( v_i(ji,jj,jl) > 0._wp ) THEN 
    423                       
     411                     ! 
    424412                     rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) 
    425413                     ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 
    426414                     ht_s  (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 
    427                       
     415                     ! 
    428416                     zdv  = v_i(ji,jj,jl) + v_s(ji,jj,jl) - zviold(ji,jj,jl) - zvsold(ji,jj,jl)   
    429                       
     417                     ! 
    430418                     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. & 
    431419                        & ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN 
    432                          
     420                        ! 
    433421                        rswitch        = MAX( 0._wp, SIGN( 1._wp, zhimax(ji,jj,jl) - epsi20 ) ) 
    434422                        a_i(ji,jj,jl)  = rswitch * ( v_i(ji,jj,jl) + v_s(ji,jj,jl) ) / MAX( zhimax(ji,jj,jl), epsi20 ) 
    435                          
     423                        ! 
    436424                        ! small correction due to *rswitch for a_i 
    437425                        v_i  (ji,jj,jl)        = rswitch * v_i  (ji,jj,jl) 
    … …  
    447435                        ENDIF 
    448436                        ! END MV MP 2016 
    449                                                  
     437                        ! 
    450438                     ENDIF 
    451                       
     439                     ! 
    452440                  ENDIF 
    453                  
     441                  ! 
    454442               END DO 
    455443            END DO 
    456444         END DO 
    457          ! ------------------------------------------------- 
    458           
    459          ! Force the upper limit of ht_i to always be < hi_max (99 m). 
    460          DO jj = 1, jpj 
     445          
     446         DO jj = 1, jpj             !--- bound ht_i to hi_max (99 m). 
    461447            DO ji = 1, jpi 
    462448               ! MV MP 2016 
    … …  
    474460            END DO 
    475461         END DO 
    476  
    477  
     462         ! 
    478463      ENDIF 
    479464          
    … …  
    482467      !------------------------------------------------------------ 
    483468      ! 
    484       at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 
    485       IF ( nn_limdyn == 1 .OR. ( ( nn_monocat == 2 ) .AND. ( jpl == 1 ) ) ) THEN ! simple conservative piling, comparable with LIM2 
     469!!gm remplace the test by, l_piling a logical compute one for all in icestp.F90 (and its declaration in ice.F90 
     470!!gm      IF ( nn_limdyn == 1 .OR. ( ( nn_monocat == 2 ) .AND. ( jpl == 1 ) ) ) THEN ! simple conservative piling, comparable with LIM2 
     471      IF( l_piling ) THEN 
     472         at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 
    486473         DO jl = 1, jpl 
    487474            DO jj = 1, jpj 
    488475               DO ji = 1, jpi 
    489                   rswitch       = MAX( 0._wp, SIGN( 1._wp, at_i(ji,jj) - epsi20 ) ) 
    490                   zda           = rswitch * MIN( rn_amax_2d(ji,jj) - at_i(ji,jj), 0._wp )  & 
    491                      &                    * a_i(ji,jj,jl) / MAX( at_i(ji,jj), epsi20 ) 
     476                  rswitch = MAX( 0._wp, SIGN( 1._wp, at_i(ji,jj) - epsi20 ) ) 
     477                  zda     = rswitch * MIN( rn_amax_2d(ji,jj) - at_i(ji,jj), 0._wp ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj), epsi20 ) 
    492478                  a_i(ji,jj,jl) = a_i(ji,jj,jl) + zda 
    493479               END DO 
    494480            END DO 
    495481         END DO 
     482!!gm better and faster coding? 
     483!         DO jl = 1, jpl 
     484!            WHERE( at_i(:,:) > epsi20 ) 
     485!               a_i(:,:,jl) = a_i(:,:,jl) * (  1._wp + MIN( rn_amax_2d(:,:) - at_i(:,:) , 0._wp ) / at_i(:,:)  ) 
     486!            END WHERE 
     487!         END DO 
     488!!gm end 
    496489      ENDIF 
    497490       
    … …  
    527520   !!   Default option         Empty Module                No sea-ice model 
    528521   !!---------------------------------------------------------------------- 
    529 CONTAINS 
    530    SUBROUTINE ice_adv        ! Empty routine 
    531    END SUBROUTINE ice_adv 
    532522#endif 
     523 
    533524   !!====================================================================== 
    534525END MODULE iceadv 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv_prather.F90

    r8409 r8486  
    55   !!====================================================================== 
    66   !! History :  LIM  ! 2008-03 (M. Vancoppenolle)  LIM-3 from LIM-2 code 
    7    !!            3.2  ! 2009-06 (F. Dupont)  correct a error in the North fold b. c. 
     7   !!            3.2  ! 2009-06 (F. Dupont)  correct a error in the North fold b.c. 
    88   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation 
    99   !!-------------------------------------------------------------------- 
    1010#if defined key_lim3 
    1111   !!---------------------------------------------------------------------- 
    12    !!   'key_lim3'                                     LIM3 sea-ice model 
    13    !!---------------------------------------------------------------------- 
    14    !!   ice_adv_x  : advection of sea ice on x axis 
    15    !!   ice_adv_y  : advection of sea ice on y axis 
    16    !!---------------------------------------------------------------------- 
    17    USE dom_oce          ! ocean domain 
    18    USE ice              ! LIM-3 variables 
     12   !!   'key_lim3'                                       LIM3 sea-ice model 
     13   !!---------------------------------------------------------------------- 
     14   !!   ice_adv_x     : advection of sea ice on x axis 
     15   !!   ice_adv_y     : advection of sea ice on y axis 
     16   !!---------------------------------------------------------------------- 
     17   USE dom_oce        ! ocean domain 
     18   USE ice            ! sea-ice variables 
    1919   ! 
    20    USE lbclnk           ! lateral boundary condition - MPP exchanges 
    21    USE in_out_manager   ! I/O manager 
    22    USE prtctl           ! Print control 
    23    USE lib_mpp          ! MPP library 
    24    USE lib_fortran      ! to use key_nosignedzero 
     20   USE lbclnk         ! lateral boundary condition - MPP exchanges 
     21   USE in_out_manager ! I/O manager 
     22   USE prtctl         ! Print control 
     23   USE lib_mpp        ! MPP library 
     24   USE lib_fortran    ! to use key_nosignedzero 
    2525 
    2626   IMPLICIT NONE 
    … …  
    3333#  include "vectopt_loop_substitute.h90" 
    3434   !!---------------------------------------------------------------------- 
    35    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     35   !! NEMO/ICE 4.0 , NEMO Consortium (2017) 
    3636   !! $Id: iceadv.F90 6746 2016-06-27 17:20:57Z clem $ 
    3737   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    4141   SUBROUTINE ice_adv_x( pdf, put , pcrh, psm , ps0 ,   & 
    4242      &                  psx, psxx, psy , psyy, psxy ) 
    43       !!--------------------------------------------------------------------- 
     43      !!---------------------------------------------------------------------- 
    4444      !!                **  routine ice_adv_x  ** 
    4545      !!   
    … …  
    5252      !! 
    5353      !! Reference:  Prather, 1986, JGR, 91, D6. 6671-6681. 
    54       !!-------------------------------------------------------------------- 
     54      !!---------------------------------------------------------------------- 
    5555      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step 
    5656      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call ice_adv_x then ice_adv_y (=1) or the opposite (=0) 
    … …  
    6868      REAL(wp), DIMENSION(jpi,jpj) ::   zfm , zfxx , zfyy  , zfxy   !  -      - 
    6969      REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q         !  -      - 
    70       !--------------------------------------------------------------------- 
     70      !----------------------------------------------------------------------- 
    7171 
    7272      ! Limitation of moments.                                            

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv_umx.F90

    r8409 r8486  
    44   !! LIM sea-ice model : sea-ice advection using the ULTIMATE-MACHO scheme 
    55   !!============================================================================== 
    6    !! History :  3.5  !  2014-11  (C. Rousset, G. Madec)  Original code 
    7    !!---------------------------------------------------------------------- 
    8  
     6   !! History :  3.6  !  2014-11  (C. Rousset, G. Madec)  Original code 
     7   !!---------------------------------------------------------------------- 
     8#if defined key_lim3 
     9   !!---------------------------------------------------------------------- 
     10   !!   'key_lim3'                                    LIM 3.0 sea-ice model 
    911   !!---------------------------------------------------------------------- 
    1012   !!   ice_adv_umx   : update the tracer trend with the 3D advection trends using a TVD scheme 
    11    !!   ultimate      : compute a tracer value at velocity points using ULTIMATE scheme at various orders 
    12    !!   macho         :  
     13   !!   ultimate_x(_y): compute a tracer value at velocity points using ULTIMATE scheme at various orders 
     14   !!   macho         : ??? 
    1315   !!   nonosc_2d     : compute monotonic tracer fluxes by a non-oscillatory algorithm  
    14    !!---------------------------------------------------------------------- 
    15 #if defined key_lim3 
    16    !!---------------------------------------------------------------------- 
    17    !!   'key_lim3' :                                  LIM 3.0 sea-ice model 
    1816   !!---------------------------------------------------------------------- 
    1917   USE phycst         ! physical constant 
    2018   USE dom_oce        ! ocean domain 
    21    USE sbc_oce, ONLY: nn_fsbc   ! ocean surface boundary condition 
    22    USE ice            ! ice variables 
     19   USE sbc_oce , ONLY : nn_fsbc   ! update frequency of surface boundary condition 
     20   USE ice            ! sea-ice variables 
    2321   ! 
    2422   USE in_out_manager ! I/O manager 
    … …  
    3331   PUBLIC   ice_adv_umx    ! routine called by iceadv.F90 
    3432       
    35    REAL(wp) ::   z1_6   = 1._wp / 6._wp   ! =1/6 
    36    REAL(wp) ::   z1_120 = 1._wp / 120._wp ! =1/120 
     33   REAL(wp) ::   z1_6   = 1._wp /   6._wp   ! =1/6 
     34   REAL(wp) ::   z1_120 = 1._wp / 120._wp   ! =1/120 
    3735 
    3836   !! * Substitutions 
    3937#  include "vectopt_loop_substitute.h90" 
    4038   !!---------------------------------------------------------------------- 
    41    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     39   !! NEMO/ICE 4.0 , NEMO Consortium (2017) 
    4240   !! $Id: iceadv_umx.F90 4499 2014-02-18 15:14:31Z timgraham $ 
    4341   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    5856      !! ** Action : - pt  the after advective tracer 
    5957      !!---------------------------------------------------------------------- 
    60       INTEGER                     , INTENT(in   )           ::   kt         ! number of iteration 
    61       REAL(wp)                    , INTENT(in   )           ::   pdt        ! tracer time-step 
    62       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   puc, pvc   ! 2 ice velocity components => u*e2 
    63       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   pubox, pvbox   ! upstream velocity 
    64       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   ptc        ! tracer content field 
     58      INTEGER                     , INTENT(in   ) ::   kt             ! number of iteration 
     59      REAL(wp)                    , INTENT(in   ) ::   pdt            ! tracer time-step 
     60      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   puc  , pvc     ! 2 ice velocity components => u*e2 
     61      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pubox, pvbox   ! upstream velocity 
     62      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ptc            ! tracer content field 
    6563      ! 
    6664      INTEGER  ::   ji, jj           ! dummy loop indices   
    … …  
    6866      REAL(wp) ::   zfp_ui, zfp_vj   !   -      - 
    6967      REAL(wp) ::   zfm_ui, zfm_vj   !   -      - 
    70       REAL(wp), DIMENSION(jpi,jpj) :: zt_ups, zfu_ups, zfv_ups, ztrd, zfu_ho, zfv_ho, zt_u, zt_v 
     68      REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ups, zfu_ho, zt_u, zt_ups 
     69      REAL(wp), DIMENSION(jpi,jpj) ::   zfv_ups, zfv_ho, zt_v, ztrd 
    7170      !!---------------------------------------------------------------------- 
    7271      ! 
    … …  
    394393      ! 
    395394      SELECT CASE (k_order ) 
    396          ! 
    397       CASE( 1 )                                                   !==  1st order central TIM  ==! (Eq. 21) 
    398          !         
     395      ! 
     396      CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21) 
    399397         DO jj = 1, jpjm1 
    400398            DO ji = 1, jpi 
    … …  
    404402         END DO 
    405403         ! 
    406       CASE( 2 )                                                   !==  2nd order central TIM  ==! (Eq. 23) 
     404      CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23) 
    407405         DO jj = 1, jpjm1 
    408406            DO ji = 1, jpi 
    … …  
    414412         CALL lbc_lnk( pt_v(:,:) , 'V',  1. ) 
    415413         ! 
    416       CASE( 3 )                                                   !==  3rd order central TIM  ==! (Eq. 24) 
    417          ! 
     414      CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24) 
    418415         DO jj = 1, jpjm1 
    419416            DO ji = 1, jpi 
    … …  
    428425         END DO 
    429426         ! 
    430       CASE( 4 )                                                   !==  4th order central TIM  ==! (Eq. 27) 
    431          ! 
     427      CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27) 
    432428         DO jj = 1, jpjm1 
    433429            DO ji = 1, jpi 
    … …  
    442438         END DO 
    443439         ! 
    444       CASE( 5 )                                                   !==  5th order central TIM  ==! (Eq. 29) 
    445          ! 
     440      CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29) 
    446441         DO jj = 1, jpjm1 
    447442            DO ji = 1, jpi 
    … …  
    569564   !!   Default option           Dummy module      NO LIM 3.0 sea-ice model 
    570565   !!---------------------------------------------------------------------- 
    571 CONTAINS 
    572    ! 
    573    SUBROUTINE ice_adv_umx     ! Dummy routine 
    574       WRITE(*,*) 'ice_adv_umx: You should not have seen this print! error?' 
    575    END SUBROUTINE ice_adv_umx 
    576566#endif 
     567 
    577568   !!====================================================================== 
    578569END MODULE iceadv_umx 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90

    r8426 r8486  
    44   !! Ocean forcing:  bulk thermohaline forcing of the ice 
    55   !!===================================================================== 
    6    !! History : 
    7    !!   NEMO     4.0  ! 2017-07  (C. Rousset) Split ice and ocean albedos 
    8    !!---------------------------------------------------------------------- 
    9    !!   ice_alb    : albedo for ice (clear and overcast skies) 
    10    !!   alb_init   : initialisation of albedo computation 
     6   !! History :  4.0  ! 2017-07  (C. Rousset) Split ice and ocean albedos 
    117   !!---------------------------------------------------------------------- 
    128#if defined key_lim3 
    139   !!---------------------------------------------------------------------- 
    14    !!   'key_lim3' :                                  LIM 3.0 sea-ice model 
    15    !!---------------------------------------------------------------------- 
    16    USE ice, ONLY : jpl 
     10   !!   'key_lim3'                                    LIM 3.0 sea-ice model 
     11   !!---------------------------------------------------------------------- 
     12   !!   ice_alb       : albedo for ice (clear and overcast skies) 
     13   !!   alb_init      : initialisation of albedo computation 
     14   !!---------------------------------------------------------------------- 
     15   USE ice     , ONLY : jpl   ! number of ice category 
    1716   USE phycst         ! physical constants 
    1817   ! 
    … …  
    2423   PRIVATE 
    2524 
    26    PUBLIC   ice_alb   ! routine called icestp.F90 
    27  
    28    REAL(wp), PUBLIC, PARAMETER ::   rn_alb_oce = 0.066   ! ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001) 
     25   PUBLIC   ice_alb   ! routine called in iceforcing.F90 and iceupdate.F90 
     26 
     27   REAL(wp), PUBLIC, PARAMETER ::   rn_alb_oce = 0.066   !: ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001) 
     28 
     29!!gm real variable name starting with a "c" NOT DOCTOR !!!! 
    2930   INTEGER  ::   albd_init = 0       ! control flag for initialization   
    30    REAL(wp) ::   c1        = 0.05    ! snow thickness (only for nn_ice_alb=0) 
    31    REAL(wp) ::   c2        = 0.10    !  "        " 
    32    REAL(wp) ::   rcloud    = 0.06    ! cloud effect on albedo (only-for nn_ice_alb=0) 
     31   REAL(wp) , PARAMETER ::   c1        = 0.05    ! snow thickness (only for nn_ice_alb=0) 
     32   REAL(wp) , PARAMETER ::   c2        = 0.10    !  "        " 
     33   REAL(wp) , PARAMETER ::   rcloud    = 0.06    ! cloud effect on albedo (only-for nn_ice_alb=0) 
     34   REAL(wp) , PARAMETER ::   r1_c1 = 1. / c1 
     35   REAL(wp) , PARAMETER ::   r1_c2 = 1. / c2 
    3336  
    34    !                             !!* namelist namsbc_alb 
     37   !                             !!* namelist namsbc_alb * 
    3538   INTEGER  ::   nn_ice_alb 
    3639   REAL(wp) ::   rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt, rn_alb_dpnd 
    3740 
    3841   !!---------------------------------------------------------------------- 
    39    !! NEMO/OPA 4.0 , NEMO Consortium (2010) 
     42   !! NEMO/ICE 4.0 , NEMO Consortium (2017) 
    4043   !! $Id: icealb.F90 8268 2017-07-03 15:01:04Z clem $ 
    4144   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    7679      !!                Brandt et al. 2005, J. Climate, vol 18 
    7780      !!                Grenfell & Perovich 2004, JGR, vol 109  
    78       !! 
    79       !!---------------------------------------------------------------------- 
    80       !! 
    81       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_ice              !  ice surface temperature (Kelvin) 
    82       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_ice              !  sea-ice thickness 
    83       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_snw              !  snow depth 
    84       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pafrac_pnd          !  melt pond relative fraction (per unit ice area) 
    85       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_pnd              !  melt pond depth 
    86       LOGICAL , INTENT(in   )                   ::   ld_pnd              !  melt ponds radiatively active or not 
    87       REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_cs           !  albedo of ice under clear    sky 
    88       REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_os           !  albedo of ice under overcast sky 
    89       ! 
    90       INTEGER  ::   ji, jj, jl                                           ! dummy loop indices 
    91       REAL(wp) ::   zswitch, z1_c1, z1_c2 
    92       REAL(wp) ::   zhref_pnd                                  
     81      !!---------------------------------------------------------------------- 
     82      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_ice       !  ice surface temperature (Kelvin) 
     83      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_ice       !  sea-ice thickness 
     84      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_snw       !  snow depth 
     85      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pafrac_pnd   !  melt pond relative fraction (per unit ice area) 
     86      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_pnd       !  melt pond depth 
     87      LOGICAL , INTENT(in   )                   ::   ld_pnd       !  melt ponds radiatively active or not 
     88      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_cs    !  albedo of ice under clear    sky 
     89      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_os    !  albedo of ice under overcast sky 
     90      ! 
     91      INTEGER  ::   ji, jj, jl                ! dummy loop indices 
     92      REAL(wp) ::   zswitch, z1_c1, z1_c2     ! local scalar 
     93      REAL(wp) ::   z1_href_pnd               !   -      -                  
    9394      REAL(wp) ::   zalb_sm, zalb_sf, zalb_st ! albedo of snow melting, freezing, total 
    9495      ! 
    95       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb, zalb_it             ! intermediate variable & albedo of ice (snow free) 
     96      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb, zalb_it   ! intermediate variable & albedo of ice (snow free) 
    9697!! MV MP 
    97       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_pnd                  ! ponded sea ice albedo 
    98       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_ice                  ! bare sea ice albedo 
    99       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_snw                  ! snow-covered sea ice albedo 
    100       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zafrac_snw                ! relative snow fraction 
    101       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zafrac_ice                ! relative ice fraction 
    102       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zafrac_pnd                ! relative ice fraction (effective) 
    103       !! 
     98      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_pnd        ! ponded sea ice albedo 
     99      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_ice        ! bare sea ice albedo 
     100      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_snw        ! snow-covered sea ice albedo 
     101      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zafrac_snw      ! relative snow fraction 
     102      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zafrac_ice      ! relative ice fraction 
     103      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zafrac_pnd      ! relative ice fraction (effective) 
    104104      !!--------------------------------------------------------------------- 
    105105 
    … …  
    120120      ENDIF 
    121121 
    122       SELECT CASE ( nn_ice_alb ) 
    123  
    124       !------------------------------------------ 
    125       !  Shine and Henderson-Sellers (1985) 
    126       !------------------------------------------ 
    127       ! NB: This parameterization is based on clear sky values 
    128  
    129       CASE( 0 ) 
    130         
    131          ! Thickness-dependent bare ice albedo 
     122      SELECT CASE ( nn_ice_alb )       ! select a parameterization 
     123      ! 
     124      !              !------------------------------------------ 
     125      CASE( 0 )      !  Shine and Henderson-Sellers (1985)     !   (based on clear sky values) 
     126         !           !------------------------------------------ 
     127         ! 
     128         !                       ! Thickness-dependent bare ice albedo 
    132129         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb 
    133130         ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = 0.472  + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 ) 
    … …  
    137134         ELSE WHERE                                       ;  zalb_it = 0.1    + 3.6    * ph_ice 
    138135         END WHERE 
    139  
    140          IF ( ld_pnd ) THEN 
    141             ! Depth-dependent ponded ice albedo 
    142             zhref_pnd = 0.05        ! Characteristic length scale for thickness dependence of ponded ice albedo, Lecomte et al (2015) 
    143             zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_it ) * EXP( - ph_pnd / zhref_pnd )  
    144  
    145             ! Snow-free ice albedo is a function of pond fraction 
    146             WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ; zalb_it = zalb_it * ( 1. - pafrac_pnd  ) + zalb_pnd * pafrac_pnd ;   END WHERE 
     136         ! 
     137         IF ( ld_pnd ) THEN      ! Depth-dependent ponded ice albedo 
     138            z1_href_pnd = 1. / 0.05    ! inverse of the characteristic length scale (Lecomte et al. 2015) 
     139            zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_it ) * EXP( - ph_pnd * z1_href_pnd )  
     140            !                          ! Snow-free ice albedo is a function of pond fraction 
     141            WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )    
     142               zalb_it = zalb_it * ( 1. - pafrac_pnd  ) + zalb_pnd * pafrac_pnd  
     143            END WHERE 
    147144         ENDIF  
    148  
     145         ! 
     146!!gm: optimization ( rn_alb_smlt - rn_alb_imlt ) * r1_c2 can be computed one for all  
     147!!gm  before the DO-loop below 
    149148         DO jl = 1, jpl 
    150149            DO jj = 1, jpj 
    151150               DO ji = 1, jpi 
    152                   ! Freezing snow 
     151                  !                    ! Freezing snow 
    153152                  ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > c2 
    154                   zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - c1 ) ) ) 
    155                   zalb_sf   = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  & 
    156                      &                           + ph_snw(ji,jj,jl) * ( rn_alb_sdry - zalb_it(ji,jj,jl) ) / c1  )   & 
    157                      &        +         zswitch   * rn_alb_sdry   
    158  
    159                   ! Melting snow 
     153                  zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - c1 ) ) ) 
     154                  zalb_sf = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  & 
     155                     &                                   + ph_snw(ji,jj,jl) * ( rn_alb_sdry - zalb_it(ji,jj,jl) ) * r1_c1  )   & 
     156                     &    +           zswitch   * rn_alb_sdry   
     157                     ! 
     158                  !                    ! Melting snow 
    160159                  ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > c2 
    161                   zswitch   = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - c2 ) ) 
    162                   zalb_sm = ( 1._wp - zswitch ) * ( rn_alb_imlt + ph_snw(ji,jj,jl) * ( rn_alb_smlt - rn_alb_imlt ) / c2 )   & 
    163                       &     +         zswitch   *   rn_alb_smlt  
     160                  zswitch = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - c2 ) ) 
     161                  zalb_sm = ( 1._wp - zswitch ) * ( rn_alb_imlt + ph_snw(ji,jj,jl) * ( rn_alb_smlt - rn_alb_imlt ) * r1_c2 )   & 
     162                     &      +         zswitch   *   rn_alb_smlt  
     163                     ! 
     164                  !                    ! Snow albedo 
     165                  zswitch =  MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )    
     166                  zalb_st =  zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf 
    164167                  ! 
    165                   ! Snow albedo 
    166                   zswitch  =  MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )    
    167                   zalb_st  =  zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf 
    168                 
    169                   ! Surface albedo 
     168                  !                    ! Surface albedo 
    170169                  zswitch             = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) ) 
    171170                  pa_ice_cs(ji,jj,jl) = zswitch * zalb_st + ( 1._wp - zswitch ) * zalb_it(ji,jj,jl) 
    … …  
    174173            END DO 
    175174         END DO 
    176  
     175         ! 
    177176         pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + rcloud       ! Oberhuber correction for overcast sky 
    178  
    179       !------------------------------------------ 
    180       !  New parameterization (2016) 
    181       !------------------------------------------ 
    182       ! NB: This parameterization is based on overcast skies values 
    183        
    184       CASE( 1 )  
    185  
    186 ! compilation of values from literature (reference overcast sky values) 
    187 !        rn_alb_sdry = 0.85      ! dry snow 
    188 !        rn_alb_smlt = 0.75      ! melting snow 
    189 !        rn_alb_idry = 0.60      ! bare frozen ice 
    190 !        rn_alb_dpnd = 0.36      ! ponded-ice overcast albedo (Lecomte et al, 2015) 
    191 !                                ! early melt pnds 0.27, late melt ponds 0.14 Grenfell & Perovich 
    192 ! Perovich et al 2002 (Sheba) => the only dataset for which all types of ice/snow were retrieved 
    193 !        rn_alb_sdry = 0.85      ! dry snow 
    194 !        rn_alb_smlt = 0.72      ! melting snow 
    195 !        rn_alb_idry = 0.65      ! bare frozen ice 
    196 ! Brandt et al 2005 (East Antarctica) 
    197 !        rn_alb_sdry = 0.87      ! dry snow 
    198 !        rn_alb_smlt = 0.82      ! melting snow 
    199 !        rn_alb_idry = 0.54      ! bare frozen ice 
    200 !  
    201          ! Computation of snow-free ice albedo  
     177         ! 
     178         ! 
     179         !           !------------------------------------------ 
     180      CASE( 1 )      !  New parameterization (2016)            !    (based on overcast sky values) 
     181         !           !------------------------------------------ 
     182         ! 
     183         !                    compilation of values from literature (reference overcast sky values) 
     184         !                          rn_alb_sdry = 0.85      ! dry snow 
     185         !                          rn_alb_smlt = 0.75      ! melting snow 
     186         !                          rn_alb_idry = 0.60      ! bare frozen ice 
     187         !                          rn_alb_dpnd = 0.36      ! ponded-ice overcast albedo (Lecomte et al, 2015) 
     188         !                                                  ! early melt pnds 0.27, late melt ponds 0.14 Grenfell & Perovich 
     189         !                    Perovich et al 2002 (Sheba) => the only dataset for which all types of ice/snow were retrieved 
     190         !                          rn_alb_sdry = 0.85      ! dry snow 
     191         !                          rn_alb_smlt = 0.72      ! melting snow 
     192         !                          rn_alb_idry = 0.65      ! bare frozen ice 
     193         !                    Brandt et al 2005 (East Antarctica) 
     194         !                          rn_alb_sdry = 0.87      ! dry snow 
     195         !                          rn_alb_smlt = 0.82      ! melting snow 
     196         !                          rn_alb_idry = 0.54      ! bare frozen ice 
     197         !  
     198         !              !--- Computation of snow-free ice albedo  
    202199         z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) )  
    203200         z1_c2 = 1. / 0.05 
    204201 
    205          ! Accounting for the ice-thickness dependency 
    206          WHERE     ( 1.5  < ph_ice                     )        ;  zalb_it = zalb 
    207          ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )        ;  zalb_it = zalb     + ( 0.18 - zalb     ) * z1_c1 *  & 
    208             &                                                                     ( LOG(1.5) - LOG(ph_ice) ) 
    209          ELSE WHERE                                             ;  zalb_it = rn_alb_oce + ( 0.18 - rn_alb_oce ) * z1_c2 * ph_ice 
    210          END WHERE 
    211  
    212          IF ( ld_pnd ) THEN 
    213             ! Depth-dependent ponded ice albedo 
    214             zhref_pnd = 0.05        ! Characteristic length scale for thickness dependence of ponded ice albedo, Lecomte et al (2015) 
    215             zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_it ) * EXP( - ph_pnd / zhref_pnd )  
    216  
    217             ! Snow-free ice albedo is weighted mean of ponded ice and bare ice contributions 
    218             WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;  zalb_it = zalb_it * ( 1. - pafrac_pnd  ) + zalb_pnd * pafrac_pnd ;  END WHERE 
     202         !              !--- Accounting for the ice-thickness dependency 
     203         WHERE     ( 1.5  < ph_ice                     )   ;   zalb_it = zalb 
     204         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )   ;   zalb_it = zalb + ( 0.18 - zalb ) * z1_c1 * ( LOG(1.5) - LOG(ph_ice) ) 
     205         ELSE WHERE                                        ;   zalb_it = rn_alb_oce + ( 0.18 - rn_alb_oce ) * z1_c2 * ph_ice 
     206         END WHERE 
     207         ! 
     208         IF ( ld_pnd ) THEN      ! Depth-dependent ponded ice albedo 
     209            z1_href_pnd = 0.05        ! inverse of the characteristic length scale (Lecomte et al. 2015) 
     210            zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_it ) * EXP( - ph_pnd * z1_href_pnd )  
     211            ! 
     212            !                    ! Snow-free ice albedo is weighted mean of ponded ice and bare ice contributions 
     213            WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice ) 
     214               zalb_it = zalb_it * ( 1. - pafrac_pnd  ) + zalb_pnd * pafrac_pnd 
     215            END WHERE 
    219216         ENDIF  
    220  
     217         ! 
     218         !              !--- Overcast sky surface albedo (accounting for snow, ice melt ponds) 
    221219         z1_c1 = 1. / 0.02 
    222220         z1_c2 = 1. / 0.03 
    223           
    224          ! Overcast sky surface albedo (accounting for snow, ice melt ponds) 
    225221         DO jl = 1, jpl 
    226222            DO jj = 1, jpj 
    … …  
    241237            END DO 
    242238         END DO 
    243  
    244          ! Clear sky surface albedo 
     239         ! 
     240         !              !--- Clear sky surface albedo 
    245241         pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 );  
    246  
    247       !--------------------------------------------------- 
    248       !  Fractional surface-based parameterization (2017) 
    249       !--------------------------------------------------- 
    250       CASE( 2 )  
    251   
    252       ! MV: I propose this formulation that is more elegant, and more easy to expand towards 
    253       !     varying pond and snow fraction. 
    254       !     Formulation 1 has issues to handle ponds and snow together that 
    255       !     can't easily be fixed. This one handles it better, I believe. 
    256  
    257           !----------------------------------------- 
    258           ! Snow, bare ice and ponded ice fractions  
    259           !----------------------------------------- 
    260           ! Specific fractions (zafrac) refer to relative area covered by snow within each ice category 
    261  
    262           !--- Effective pond fraction (for now, we prevent melt ponds and snow at the same time) 
    263           zafrac_pnd = 0._wp 
    264           IF ( ld_pnd ) THEN   
    265              WHERE( ph_snw == 0._wp ) ;  zafrac_pnd = pafrac_pnd ;  END WHERE  ! Snow fully "shades" melt ponds 
    266           ENDIF          
    267  
    268           !--- Specific snow fraction (for now, prescribed) 
    269           WHERE     ( ph_snw > 0._wp     ) ;  zafrac_snw = 1. 
    270           ELSE WHERE                       ;  zafrac_snw = 0. 
    271           END WHERE 
    272   
    273           !--- Specific ice fraction  
    274           zafrac_ice = 1. - zafrac_snw - zafrac_pnd 
    275   
    276           !-------------------------------------------------- 
    277           ! Snow-covered, pond-covered, and bare ice albedos 
    278           !-------------------------------------------------- 
    279           ! Bare ice albedo 
    280           z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) )  
    281           z1_c2 = 1. / 0.05 
    282           WHERE     ( 1.5  < ph_ice                     )  ;  zalb_ice = zalb 
    283           ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_ice = zalb     + ( 0.18 - zalb     ) * z1_c1 *  & 
     242         ! 
     243         ! 
     244         !           !------------------------------------------ 
     245      CASE( 2 )      !  Fractional surface-based parameterization (2017) 
     246         !           !------------------------------------------ 
     247         !              MV: I propose this formulation that is more elegant, and more easy to expand towards 
     248         !              varying pond and snow fraction. 
     249         !              Formulation 1 has issues to handle ponds and snow together that can't easily be fixed.  
     250         !              This one handles it better, I believe. 
     251         ! 
     252         !                 !==  Snow, bare ice and ponded ice fractions  ==! 
     253         ! 
     254         !                       ! Specific fractions (zafrac) refer to relative area covered by snow within each ice category 
     255         ! 
     256         !                       !--- Effective pond fraction (for now, we prevent melt ponds and snow at the same time) 
     257         zafrac_pnd = 0._wp 
     258         IF ( ld_pnd ) THEN   
     259            WHERE( ph_snw == 0._wp ) ;  zafrac_pnd = pafrac_pnd ;  END WHERE  ! Snow fully "shades" melt ponds 
     260         ENDIF          
     261         ! 
     262         !                       !--- Specific snow fraction (for now, prescribed) 
     263         WHERE     ( ph_snw > 0._wp     ) ;  zafrac_snw = 1. 
     264         ELSE WHERE                       ;  zafrac_snw = 0. 
     265         END WHERE 
     266         ! 
     267         !                       !--- Specific ice fraction  
     268         zafrac_ice = 1. - zafrac_snw - zafrac_pnd 
     269         ! 
     270         !                 !==  Snow-covered, pond-covered, and bare ice albedos  ==! 
     271         ! 
     272         !                       !--- Bare ice albedo 
     273         z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) )  
     274         z1_c2 = 1. / 0.05 
     275         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_ice = zalb 
     276         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_ice = zalb     + ( 0.18 - zalb     ) * z1_c1 *  & 
    284277            &                                                                       ( LOG(1.5) - LOG(ph_ice) ) 
    285           ELSE WHERE                                       ;  zalb_ice = rn_alb_oce + ( 0.18 - rn_alb_oce ) * z1_c2 * ph_ice 
    286           END WHERE 
    287  
    288           ! Snow-covered ice albedo (freezing, melting cases) 
    289           z1_c1 = 1. / 0.02 
    290           z1_c2 = 1. / 0.03 
    291           
    292           WHERE( pt_ice < rt0_snow ) ; zalb_snw = rn_alb_sdry - ( rn_alb_sdry - zalb_ice ) * EXP( - ph_snw * z1_c1 ); 
    293           ELSE WHERE                 ; zalb_snw = rn_alb_smlt - ( rn_alb_smlt - zalb_ice ) * EXP( - ph_snw * z1_c2 ); 
    294           END WHERE 
    295  
    296           ! Depth-dependent ponded ice albedo 
    297           IF ( ld_pnd ) THEN 
    298              zhref_pnd = 0.05        ! Characteristic length scale for thickness dependence of ponded ice albedo, Lecomte et al (2015) 
    299              zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_ice ) * EXP( - ph_pnd / zhref_pnd )  
    300           ELSE 
    301              zalb_pnd  = rn_alb_dpnd 
    302           ENDIF 
    303  
    304           ! Surface albedo is weighted mean of snow, ponds and bare ice contributions 
    305           pa_ice_os = zafrac_snw * zalb_snw  +  zafrac_pnd * zalb_pnd  +  zafrac_ice * zalb_ice 
    306            
    307           pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 ) 
     278         ELSE WHERE                                       ;  zalb_ice = rn_alb_oce + ( 0.18 - rn_alb_oce ) * z1_c2 * ph_ice 
     279         END WHERE 
     280         ! 
     281         z1_c1 = 1. / 0.02       !--- Snow-covered ice albedo (freezing, melting cases) 
     282         z1_c2 = 1. / 0.03 
     283         ! 
     284         WHERE( pt_ice < rt0_snow ) ; zalb_snw = rn_alb_sdry - ( rn_alb_sdry - zalb_ice ) * EXP( - ph_snw * z1_c1 ); 
     285         ELSE WHERE                 ; zalb_snw = rn_alb_smlt - ( rn_alb_smlt - zalb_ice ) * EXP( - ph_snw * z1_c2 ); 
     286         END WHERE 
     287         ! 
     288         IF ( ld_pnd ) THEN      !--- Depth-dependent ponded ice albedo 
     289            z1_href_pnd = 0.05        ! inverse of the characteristic length scale (Lecomte et al. 2015) 
     290            zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_ice ) * EXP( - ph_pnd * z1_href_pnd )  
     291         ELSE 
     292            zalb_pnd  = rn_alb_dpnd 
     293         ENDIF 
     294         !                       !--- Surface albedo is weighted mean of snow, ponds and bare ice contributions 
     295         pa_ice_os = zafrac_snw * zalb_snw  +  zafrac_pnd * zalb_pnd  +  zafrac_ice * zalb_ice 
     296         pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 ) 
    308297 
    309298      END SELECT 
    … …  
    311300   END SUBROUTINE ice_alb 
    312301 
     302 
    313303   SUBROUTINE alb_init 
    314304      !!---------------------------------------------------------------------- 
    … …  
    320310      !!---------------------------------------------------------------------- 
    321311      INTEGER  ::   ios                 ! Local integer output status for namelist read 
     312      !! 
    322313      NAMELIST/namicealb/ nn_ice_alb, rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt, rn_alb_dpnd 
    323314      !!---------------------------------------------------------------------- 
    … …  
    353344   !!   Default option           Dummy module      NO LIM 3.0 sea-ice model 
    354345   !!---------------------------------------------------------------------- 
    355 CONTAINS 
    356    ! 
    357    SUBROUTINE ice_alb     ! Dummy routine 
    358       WRITE(*,*) 'ice_alb: You should not have seen this print! error?' 
    359    END SUBROUTINE ice_alb 
    360346#endif 
    361347 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icecor.F90

    r8426 r8486  
    99#if defined key_lim3 
    1010   !!---------------------------------------------------------------------- 
    11    !!   'key_lim3'                                      LIM3 sea-ice model 
    12    !!---------------------------------------------------------------------- 
    13    !!    ice_cor   : computes update of sea-ice global variables from trend terms 
    14    !!---------------------------------------------------------------------- 
    15    USE dom_oce 
    16    USE phycst          ! physical constants 
    17    USE ice 
    18    USE ice1D           ! LIM thermodynamic sea-ice variables 
    19    USE iceitd 
    20    USE icevar 
    21    USE icectl          ! control prints 
     11   !!   'key_lim3'                                       LIM3 sea-ice model 
     12   !!---------------------------------------------------------------------- 
     13   !!    ice_cor      : computes update of sea-ice global variables from trend terms 
     14   !!---------------------------------------------------------------------- 
     15   USE dom_oce        ! ocean domain 
     16   USE phycst         ! physical constants 
     17   USE ice            ! sea-ice: variable 
     18   USE ice1D          ! sea-ice: thermodynamic sea-ice variables 
     19   USE iceitd         ! sea-ice: rebining 
     20   USE icevar         ! sea-ice: operations 
     21   USE icectl         ! sea-ice: control prints 
    2222   ! 
    23    USE in_out_manager  ! I/O manager 
    24    USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    25    USE lbclnk          ! lateral boundary condition - MPP link 
    26    USE lib_mpp         ! MPP library 
    27    USE timing          ! Timing 
     23   USE in_out_manager ! I/O manager 
     24   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     25   USE lbclnk         ! lateral boundary condition - MPP link 
     26   USE lib_mpp        ! MPP library 
     27   USE timing         ! Timing 
    2828 
    2929   IMPLICIT NONE 
    3030   PRIVATE 
    3131 
    32    PUBLIC   ice_cor 
     32   PUBLIC   ice_cor   ! called by icestp.F90 
    3333 
    3434   !! * Substitutions 
    3535#  include "vectopt_loop_substitute.h90" 
    3636   !!---------------------------------------------------------------------- 
    37    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     37   !! NEMO/ICE 4.0 , NEMO Consortium (2017) 
    3838   !! $Id: icecor.F90 8378 2017-07-26 13:55:59Z clem $ 
    3939   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    4141CONTAINS 
    4242 
    43    SUBROUTINE ice_cor( kt , kn ) 
    44       !!------------------------------------------------------------------- 
     43   SUBROUTINE ice_cor( kt, kn ) 
     44      !!---------------------------------------------------------------------- 
    4545      !!               ***  ROUTINE ice_cor  *** 
    4646      !!                
    47       !! ** Purpose :  Computes update of sea-ice global variables at  
    48       !!               the end of the dynamics. 
    49       !!                 
    50       !!--------------------------------------------------------------------- 
     47      !! ** Purpose :   Computes corrections on sea-ice global variables at  
     48      !!              the end of the dynamics. 
     49      !!---------------------------------------------------------------------- 
    5150      INTEGER, INTENT(in) ::   kt    ! number of iteration 
    5251      INTEGER, INTENT(in) ::   kn    ! 1 = after dyn ; 2 = after thermo 
     52      ! 
    5353      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    54       REAL(wp) ::   zsal 
    55       REAL(wp) ::   zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    56       !!------------------------------------------------------------------- 
    57       IF( nn_timing == 1 )  CALL timing_start('icecor') 
    58  
     54      REAL(wp) ::   zsal, zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b, zzc 
     55      !!---------------------------------------------------------------------- 
     56      IF( nn_timing == 1 )   CALL timing_start('icecor') 
     57      ! 
    5958      IF( kt == nit000 .AND. lwp .AND. kn == 2 ) THEN 
    6059         WRITE(numout,*) 
    61          WRITE(numout,*)' icecor '  
    62          WRITE(numout,*)' ~~~~~~ ' 
    63       ENDIF 
    64  
    65       ! conservation test 
    66       IF( ln_limdiachk ) CALL ice_cons_hsm(0, 'icecor', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    67  
    68       !---------------------------------------------------------------------- 
    69       ! Constrain the thickness of the smallest category above himin 
    70       !---------------------------------------------------------------------- 
    71       IF( kn == 2 ) THEN 
    72           
     60         WRITE(numout,*)' icecor :  correct sea ice variables if out of bounds '  
     61         WRITE(numout,*)' ~~~~~~~' 
     62      ENDIF 
     63      !                             !--- conservation test 
     64      IF( ln_limdiachk )   CALL ice_cons_hsm(0, 'icecor', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     65      ! 
     66      ! 
     67      !                             !----------------------------------------------------- 
     68      IF( kn == 2 ) THEN            !  thickness of the smallest category above himin    ! 
     69         !                          !----------------------------------------------------- 
     70         ! 
    7371         DO jj = 1, jpj  
    7472            DO ji = 1, jpi 
     73!!gm replace this 
    7574               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,1) - epsi20 ) )   !0 if no ice and 1 if yes 
    7675               ht_i(ji,jj,1) = v_i (ji,jj,1) / MAX( a_i(ji,jj,1) , epsi20 ) * rswitch 
     76!!gm by more readable coding (not slower coding since already a IF in the loop): 
     77!               IF( a_i(ji,jj,1) >= epsi20 )   ht_i(ji,jj,1) = v_i (ji,jj,1) / a_i(ji,jj,1) 
     78!!gm 
    7779               IF( v_i(ji,jj,1) > 0._wp .AND. ht_i(ji,jj,1) < rn_himin ) THEN 
    7880                  a_i (ji,jj,1) = a_i (ji,jj,1) * ht_i(ji,jj,1) / rn_himin 
    … …  
    8082            END DO 
    8183         END DO 
    82  
    83       ENDIF 
    84           
    85       !---------------------------------------------------- 
    86       ! ice concentration should not exceed amax  
    87       !----------------------------------------------------- 
    88       at_i(:,:) = 0._wp 
    89       DO jl = 1, jpl 
     84         ! 
     85      ENDIF 
     86 
     87      !                             !----------------------------------------------------- 
     88      at_i(:,:) = a_i(:,:,1)        !  ice concentration should not exceed amax          ! 
     89      DO jl = 2, jpl                !----------------------------------------------------- 
    9090         at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    9191      END DO 
    92  
     92      ! 
     93!!gm   Question   it seams to me that we have the following equality (dropping the "(ji,jj)": 
     94!      ( 1. - ( 1. - rn_amax_2d / at_i ) ) =  ( 1. - ( at_i - rn_amax_2d ) / at_i ) 
     95!                                          =  ( at_i - ( at_i - rn_amax_2d ) ) / at_i 
     96!                                          =  ( + rn_amax_2d  ) / at_i 
     97!                                          =  rn_amax_2d / at_i 
     98!     No ?  if yes see "!!gm   better" juste below  
     99!gm 
    93100      DO jl  = 1, jpl 
    94101         DO jj = 1, jpj 
    95102            DO ji = 1, jpi 
    96103               IF( at_i(ji,jj) > rn_amax_2d(ji,jj) .AND. a_i(ji,jj,jl) > 0._wp ) THEN 
    97                   a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax_2d(ji,jj) / at_i(ji,jj) ) ) 
     104                  a_i(ji,jj,jl) = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax_2d(ji,jj) / at_i(ji,jj) ) ) 
     105!!gm   better:    a_i(ji,jj,jl) = a_i(ji,jj,jl) * rn_amax_2d(ji,jj) / at_i(ji,jj)            
    98106               ENDIF 
    99107            END DO 
    100108         END DO 
    101109      END DO 
     110!!gm  Other question:  why testing a_i(ji,jj,jl) > 0._wp ?   a_i is >=0, a multiplication by 0 does not change the results.... 
     111!!gm  so at the end, the loop can be recoded without IF as: 
     112!      WHERE( at_i(:,:) > rn_amax_2d(:,:) ) 
     113!         DO jl  = 1, jpl 
     114!            a_i(:,:,jl) = a_i(:,:,jl) * MAX( rn_amax_2d(:,:), at_i(:,:) ) / at_i(:,:) 
     115!         END DO 
     116!      END WHERE 
     117!!gm  No? 
    102118     
    103       !--------------------- 
    104       ! Ice salinity bounds 
    105       !--------------------- 
    106       IF (  nn_icesal == 2  ) THEN  
    107          DO jl = 1, jpl 
     119      !                             !----------------------------------------------------- 
     120      IF (  nn_icesal == 2  ) THEN  !  Ice salinity bounds                               ! 
     121      !                             !----------------------------------------------------- 
     122         zzc = rhoic * r1_rdtice 
     123         DO jl = 1, jpl                  ! salinity stays in bounds [Simin,Simax] 
    108124            DO jj = 1, jpj  
    109125               DO ji = 1, jpi 
    110                   zsal            = smv_i(ji,jj,jl) 
    111                   ! salinity stays in bounds 
    112                   rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) ) 
    113                   smv_i(ji,jj,jl) = rswitch * MAX( MIN( rn_simax * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), rn_simin * v_i(ji,jj,jl) ) 
     126                  zsal = smv_i(ji,jj,jl) 
     127                  smv_i(ji,jj,jl) = MIN(  MAX( rn_simin*v_i(ji,jj,jl) , smv_i(ji,jj,jl) ) , rn_simax*v_i(ji,jj,jl)  ) 
    114128                  ! associated salt flux 
    115                   sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 
     129                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * zzc 
    116130               END DO 
    117131            END DO 
    … …  
    119133      ENDIF 
    120134 
    121       !---------------------------------------------------- 
    122       ! Rebin categories with thickness out of bounds 
    123       !---------------------------------------------------- 
    124       IF ( jpl > 1 ) CALL ice_itd_reb 
    125  
    126       !----------------- 
    127       ! zap small values 
    128       !----------------- 
    129       CALL ice_var_zapsmall 
    130  
    131       !---------------------------------------------- 
    132       ! Ice drift. Corrections to avoid wrong values 
    133       !---------------------------------------------- 
    134       IF( kn == 2 ) THEN 
    135          DO jj = 2, jpjm1 
     135      !                             !----------------------------------------------------- 
     136      !                             !  Rebin categories with thickness out of bounds     ! 
     137      !                             !----------------------------------------------------- 
     138      IF ( jpl > 1 )   CALL ice_itd_reb 
     139 
     140      !                             !----------------------------------------------------- 
     141      CALL ice_var_zapsmall         !  Zap small values                                  ! 
     142      !                             !----------------------------------------------------- 
     143 
     144      !                             !----------------------------------------------------- 
     145      IF( kn == 2 ) THEN            !  Ice drift case: Corrections to avoid wrong values ! 
     146         DO jj = 2, jpjm1           !----------------------------------------------------- 
    136147            DO ji = 2, jpim1 
    137                IF ( at_i(ji,jj) == 0._wp ) THEN ! what to do if there is no ice 
    138                   IF ( at_i(ji+1,jj) == 0._wp ) u_ice(ji,jj)   = 0._wp ! right side 
    139                   IF ( at_i(ji-1,jj) == 0._wp ) u_ice(ji-1,jj) = 0._wp ! left side 
    140                   IF ( at_i(ji,jj+1) == 0._wp ) v_ice(ji,jj)   = 0._wp ! upper side 
    141                   IF ( at_i(ji,jj-1) == 0._wp ) v_ice(ji,jj-1) = 0._wp ! bottom side 
     148               IF ( at_i(ji,jj) == 0._wp ) THEN    ! what to do if there is no ice 
     149                  IF ( at_i(ji+1,jj) == 0._wp )   u_ice(ji  ,jj) = 0._wp  ! right side 
     150                  IF ( at_i(ji-1,jj) == 0._wp )   u_ice(ji-1,jj) = 0._wp  ! left side 
     151                  IF ( at_i(ji,jj+1) == 0._wp )   v_ice(ji,jj  ) = 0._wp  ! upper side 
     152                  IF ( at_i(ji,jj-1) == 0._wp )   v_ice(ji,jj-1) = 0._wp  ! bottom side 
    142153               ENDIF 
    143154            END DO 
    144155         END DO 
    145          !lateral boundary conditions 
    146          CALL lbc_lnk_multi( u_ice, 'U', -1., v_ice, 'V', -1. ) 
    147          !mask velocities 
    148          u_ice(:,:) = u_ice(:,:) * umask(:,:,1) 
     156         CALL lbc_lnk_multi( u_ice, 'U', -1., v_ice, 'V', -1. )            ! lateral boundary conditions 
     157         ! 
     158!!gm  I think masking here is unnecessary, u_ice already masked and we only introduce zeros in the field 
     159         u_ice(:,:) = u_ice(:,:) * umask(:,:,1)                            ! mask velocities 
    149160         v_ice(:,:) = v_ice(:,:) * vmask(:,:,1) 
    150161      ENDIF 
    151        
    152       ! ------------------------------------------------- 
    153       ! Diagnostics 
    154       ! ------------------------------------------------- 
    155       IF( kn == 1 ) THEN 
     162 
     163!!gm I guess the trends are only out on demand  
     164!!   So please, only do this is it exite an iom_use of on a these variables 
     165!!   furthermore, only allocate the diag_ arrays in this case  
     166!!   and do the iom_put here so that it is only a local allocation 
     167!!gm  
     168      !                             !----------------------------------------------------- 
     169      SELECT CASE( kn )             !  Diagnostics                                       ! 
     170      !                             !----------------------------------------------------- 
     171      CASE( 1 )                        !--- dyn trend diagnostics 
    156172         DO jl  = 1, jpl 
    157173            afx_dyn(:,:) = afx_dyn(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice 
    158174         END DO 
    159           
     175         ! 
     176!!gm   here I think the number of ice cat is too small to use a SUM instruction... 
    160177         DO jj = 1, jpj 
    161178            DO ji = 1, jpi             
    162                ! heat content variation (W.m-2) 
    163                diag_heat(ji,jj) = - ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  &  
    164                   &                   SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    & 
    165                   &                 ) * r1_rdtice 
    166                ! salt, volume 
     179               !                 ! heat content variation (W.m-2) 
     180               diag_heat(ji,jj) = - (  SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) )    &  
     181                  &                  + SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )  ) * r1_rdtice 
     182               !                 ! salt, volume 
    167183               diag_smvi(ji,jj) = SUM( smv_i(ji,jj,:) - smv_i_b(ji,jj,:) ) * rhoic * r1_rdtice 
    168184               diag_vice(ji,jj) = SUM( v_i  (ji,jj,:) - v_i_b  (ji,jj,:) ) * rhoic * r1_rdtice 
    … …  
    170186            END DO 
    171187         END DO 
    172           
    173       ELSEIF( kn == 2 ) THEN 
    174           
     188         ! 
     189      CASE( 2 )                        !--- thermo trend diagnostics & ice aging 
     190         ! 
    175191         DO jl  = 1, jpl 
    176             oa_i(:,:,jl) = oa_i(:,:,jl) + a_i(:,:,jl) * rdt_ice          ! ice natural aging 
    177             afx_thd(:,:) = afx_thd(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice 
    178          END DO 
    179          afx_tot = afx_thd + afx_dyn 
    180           
     192            oa_i(:,:,jl) = oa_i(:,:,jl) +   a_i(:,:,jl)                   * rdt_ice       ! ice natural aging incrementation 
     193            afx_thd(:,:) = afx_thd(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice     ! thermo tendancy 
     194         END DO 
     195         afx_tot(:,:) = afx_thd(:,:) + afx_dyn(:,:) 
     196         ! 
     197!!gm   here I think the number of ice cat is too small to use a SUM instruction... 
    181198         DO jj = 1, jpj 
    182199            DO ji = 1, jpi             
    183                ! heat content variation (W.m-2) 
    184                diag_heat(ji,jj) = diag_heat(ji,jj) -  & 
    185                   &               ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  &  
    186                   &                 SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    & 
    187                   &               ) * r1_rdtice    
    188                ! salt, volume 
     200               !                 ! heat content variation (W.m-2) 
     201               diag_heat(ji,jj) = diag_heat(ji,jj) - (  SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) )    &  
     202                  &                                   + SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )  ) * r1_rdtice 
     203               !                 ! salt, volume 
    189204               diag_smvi(ji,jj) = diag_smvi(ji,jj) + SUM( smv_i(ji,jj,:) - smv_i_b(ji,jj,:) ) * rhoic * r1_rdtice 
    190205               diag_vice(ji,jj) = diag_vice(ji,jj) + SUM( v_i  (ji,jj,:) - v_i_b  (ji,jj,:) ) * rhoic * r1_rdtice 
    … …  
    192207            END DO 
    193208         END DO 
    194           
    195       ENDIF 
    196        
    197       ! conservation test 
    198       IF( ln_limdiachk ) CALL ice_cons_hsm(1, 'icecor', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    199  
    200       ! control prints 
    201       IF( ln_ctl )       CALL ice_prt3D( 'icecor' ) 
    202       IF( ln_limctl .AND. kn == 2 )  & 
    203          &               CALL ice_prt( kt, iiceprt, jiceprt, 2, ' - Final state - ' ) 
    204     
    205       IF( nn_timing == 1 )  CALL timing_stop('icecor') 
    206  
     209         ! 
     210      END SELECT 
     211      ! 
     212      !                                !--- conservation test 
     213      IF( ln_limdiachk )   CALL ice_cons_hsm(1, 'icecor', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     214      ! 
     215      !                                !--- control prints 
     216      IF( ln_ctl )                    CALL ice_prt3D( 'icecor' ) 
     217      IF( ln_limctl .AND. kn == 2 )   CALL ice_prt( kt, iiceprt, jiceprt, 2, ' - Final state - ' ) 
     218      ! 
     219      IF( nn_timing == 1 )   CALL timing_stop('icecor') 
     220      ! 
    207221   END SUBROUTINE ice_cor 
    208222 
     223#else 
     224   !!---------------------------------------------------------------------- 
     225   !!   Default option           Dummy module      NO LIM 3.0 sea-ice model 
     226   !!---------------------------------------------------------------------- 
    209227#endif 
    210228 
     229   !!====================================================================== 
    211230END MODULE icecor 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90

    r8426 r8486  
    99#if defined key_lim3 
    1010   !!---------------------------------------------------------------------- 
    11    !!   'key_lim3'                                      LIM3 sea-ice model 
     11   !!   'key_lim3'                                       LIM3 sea-ice model 
    1212   !!---------------------------------------------------------------------- 
    1313   !!    ice_ctl   : control prints in case of crash 
    … …  
    1515   !!    ice_prt3D : control prints of ice arrays 
    1616   !!---------------------------------------------------------------------- 
    17    USE oce             ! ocean dynamics and tracers 
    18    USE dom_oce         ! ocean space and time domain 
    19    USE ice             ! LIM-3: ice variables 
    20    USE ice1D           ! LIM-3: thermodynamical variables 
    21    USE sbc_oce         ! Surface boundary condition: ocean fields 
    22    USE sbc_ice         ! Surface boundary condition: ice   fields 
    23    USE phycst          ! Define parameters for the routines 
     17   USE oce            ! ocean dynamics and tracers 
     18   USE dom_oce        ! ocean space and time domain 
     19   USE ice            ! LIM-3: ice variables 
     20   USE ice1D          ! LIM-3: thermodynamical variables 
     21   USE sbc_oce        ! Surface boundary condition: ocean fields 
     22   USE sbc_ice        ! Surface boundary condition: ice   fields 
     23   USE phycst         ! Define parameters for the routines 
    2424   ! 
    25    USE lib_mpp         ! MPP library 
    26    USE timing          ! Timing 
    27    USE in_out_manager  ! I/O manager 
    28    USE prtctl          ! Print control 
    29    USE lib_fortran     !  
     25   USE lib_mpp        ! MPP library 
     26   USE timing         ! Timing 
     27   USE in_out_manager ! I/O manager 
     28   USE prtctl         ! Print control 
     29   USE lib_fortran    !  
    3030 
    3131   IMPLICIT NONE 
    … …  
    4141#  include "vectopt_loop_substitute.h90" 
    4242   !!---------------------------------------------------------------------- 
    43    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     43   !! NEMO/ICE 4.0 , NEMO Consortium (2017) 
    4444   !! $Id: icectl.F90 5043 2015-01-28 16:44:18Z clem $ 
    4545   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4646   !!---------------------------------------------------------------------- 
    47  
    4847CONTAINS 
    4948 
    5049   SUBROUTINE ice_cons_hsm( icount, cd_routine, zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b ) 
    51       !!-------------------------------------------------------------------------------------------------------- 
    52       !!                                        ***  ROUTINE ice_cons_hsm *** 
     50      !!---------------------------------------------------------------------- 
     51      !!                       ***  ROUTINE ice_cons_hsm *** 
    5352      !! 
    5453      !! ** Purpose : Test the conservation of heat, salt and mass for each ice routine 
    … …  
    6160      !!              For salt and heat thresholds, ice is considered to have a salinity of 10  
    6261      !!              and a heat content of 3e5 J/kg (=latent heat of fusion)  
    63       !!-------------------------------------------------------------------------------------------------------- 
    64       INTEGER         , INTENT(in)    :: icount        ! determine wether this is the beggining of the routine (0) or the end (1) 
    65       CHARACTER(len=*), INTENT(in)    :: cd_routine    ! name of the routine 
    66       REAL(wp)        , INTENT(inout) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
     62      !!---------------------------------------------------------------------- 
     63      INTEGER         , INTENT(in)    ::   icount        ! called at: =0 the begining of the routine, =1  the end 
     64      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      !! 
    6767      REAL(wp)                        :: zvi,   zsmv,   zei,   zfs,   zfw,   zft 
    6868      REAL(wp)                        :: zvmin, zamin, zamax  
    … …  
    7070      REAL(wp)                        :: zarea, zv_sill, zs_sill, zh_sill 
    7171      REAL(wp), PARAMETER             :: zconv = 1.e-9 ! convert W to GW and kg to Mt 
    72  
     72      !!---------------------------------------------------------------------- 
     73 
     74!!gm  Note that glo_sum for a 2D or 3D array use a multiplication by tmask_i(ji,jj) 
     75!!    so below  the  * tmask(:,:,1) is useless   ===>> I have removed them 
     76!!    I also move the conversion factor from then glo_sum argument (become a single multiplication 
     77  
    7378      IF( icount == 0 ) THEN 
    74  
    75          ! salt flux 
     79         !                          ! salt flux 
    7680         zfs_b  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
    7781            &                  sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)    & 
    78             &                ) *  e1e2t(:,:) * tmask(:,:,1) * zconv ) 
    79  
    80          ! water flux 
     82            &                ) *  e1e2t(:,:) ) * zconv  
     83         ! 
     84         !                          ! water flux 
    8185         zfw_b  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:)     + wfx_sum(:,:)     + wfx_sni(:,:)     +                     & 
    8286            &                  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(:,:) * tmask(:,:,1) * zconv ) 
    85  
    86          ! heat flux 
     87            &                  wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)    & 
     88            &                ) * e1e2t(:,:) ) * zconv 
     89         ! 
     90         !                          ! heat flux 
    8791         zft_b  = glob_sum(  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  &  
    8892            &                - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   & 
    89             &                ) *  e1e2t(:,:) * tmask(:,:,1) * zconv ) 
    90  
    91          zvi_b  = glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t * tmask(:,:,1) * zconv ) 
    92  
    93          zsmv_b = glob_sum( SUM( smv_i * rhoic            , dim=3 ) * e1e2t * tmask(:,:,1) * 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 )    & 
    97                             ) * e1e2t * tmask(:,:,1) * zconv ) 
     93            &                ) *  e1e2t(:,:) ) * zconv 
     94 
     95         zvi_b  = glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t * zconv ) 
     96 
     97         zsmv_b = glob_sum( SUM( smv_i * rhoic            , dim=3 ) * e1e2t * zconv ) 
     98 
     99         zei_b  = glob_sum( (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )     & 
     100            &                + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )   ) * e1e2t ) * zconv 
    98101 
    99102      ELSEIF( icount == 1 ) THEN 
    … …  
    102105         zfs  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
    103106            &                sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)    &  
    104             &              ) * e1e2t(:,:) * tmask(:,:,1) * zconv ) - zfs_b 
     107            &              ) * e1e2t(:,:) ) * zconv - zfs_b 
    105108 
    106109         ! water flux 
    … …  
    108111            &                wfx_opw(:,:) + wfx_res(:,:)     + wfx_dyn(:,:)     + wfx_lam(:,:)     + wfx_ice_sub(:,:) +  & 
    109112            &                wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)        & 
    110             &              ) * e1e2t(:,:) * tmask(:,:,1) * zconv ) - zfw_b 
     113            &              ) * e1e2t(:,:) ) * zconv - zfw_b 
    111114 
    112115         ! heat flux 
    113116         zft  = glob_sum(  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  &  
    114117            &              - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   & 
    115             &              ) * e1e2t(:,:) * tmask(:,:,1) * zconv ) - zft_b 
     118            &              ) * e1e2t(:,:) ) * zconv - zft_b 
    116119  
    117120         ! outputs 
    118          zvi  = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) & 
    119             &                    * e1e2t * tmask(:,:,1) * zconv ) - zvi_b ) * r1_rdtice - zfw ) * rday 
    120  
    121          zsmv = ( ( glob_sum( SUM( smv_i * rhoic            , dim=3 ) & 
    122             &                    * e1e2t * tmask(:,:,1) * zconv ) - zsmv_b ) * r1_rdtice + zfs ) * rday 
    123  
    124          zei  = ( glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) +  & 
    125             &                 SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )    & 
    126             &                ) * e1e2t * tmask(:,:,1) * zconv ) - zei_b ) * r1_rdtice + zft 
     121         zvi  = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t  ) * zconv & 
     122            &       - zvi_b ) * r1_rdtice - zfw ) * rday 
     123 
     124         zsmv = ( ( glob_sum( SUM( smv_i * rhoic            , dim=3 ) * e1e2t ) * zconv & 
     125            &       - zsmv_b ) * r1_rdtice + zfs ) * rday 
     126 
     127         zei  = ( glob_sum( (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )   & 
     128            &                + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv   & 
     129            &   - zei_b ) * r1_rdtice + zft 
    127130 
    128131         ! zvtrp and zetrp must be close to 0 if the advection scheme is conservative 
    129          zvtrp = glob_sum( ( diag_trp_vi * rhoic + diag_trp_vs * rhosn ) * e1e2t * tmask(:,:,1) * zconv ) * rday  
    130          zetrp = glob_sum( ( diag_trp_ei         + diag_trp_es         ) * e1e2t * tmask(:,:,1) * zconv ) 
     132         zvtrp = glob_sum( ( diag_trp_vi * rhoic + diag_trp_vs * rhosn ) * e1e2t  ) * zconv * rday  
     133         zetrp = glob_sum( ( diag_trp_ei         + diag_trp_es         ) * e1e2t  ) * zconv 
    131134 
    132135         zvmin = glob_min( v_i ) 
    … …  
    135138 
    136139         ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice)  
    137          zarea   = glob_sum( SUM( a_i + epsi10, dim=3 ) * e1e2t * zconv ) ! in 1.e9 m2 
     140         zarea   = glob_sum( SUM( a_i + epsi10, dim=3 ) * e1e2t ) * zconv ! in 1.e9 m2 
    138141         zv_sill = zarea * 2.5e-5 
    139142         zs_sill = zarea * 25.e-5 
    … …  
    156159            IF (      zamin  < -epsi10 ) WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin 
    157160         ENDIF 
    158  
     161         ! 
    159162      ENDIF 
    160163 
    … …  
    163166 
    164167   SUBROUTINE ice_cons_final( cd_routine ) 
    165       !!--------------------------------------------------------------------------------------------------------- 
    166       !!                                   ***  ROUTINE ice_cons_final *** 
     168      !!---------------------------------------------------------------------- 
     169      !!                     ***  ROUTINE ice_cons_final *** 
    167170      !! 
    168171      !! ** Purpose : Test the conservation of heat, salt and mass at the end of each ice time-step 
    … …  
    174177      !!              For salt and heat thresholds, ice is considered to have a salinity of 10  
    175178      !!              and a heat content of 3e5 J/kg (=latent heat of fusion)  
    176       !!-------------------------------------------------------------------------------------------------------- 
     179      !!---------------------------------------------------------------------- 
    177180      CHARACTER(len=*), INTENT(in)    :: cd_routine    ! name of the routine 
    178181      REAL(wp)                        :: zhfx, zsfx, zvfx 
    179182      REAL(wp)                        :: zarea, zv_sill, zs_sill, zh_sill 
    180183      REAL(wp), PARAMETER             :: zconv = 1.e-9 ! convert W to GW and kg to Mt 
     184      !!---------------------------------------------------------------------- 
    181185 
    182186      ! heat flux 
    183187      zhfx  = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es   & 
    184188      !  &              - SUM( qevap_ice * a_i_b, dim=3 )                           & !!clem: I think this line must be commented (but need check) 
    185          &              ) * e1e2t * tmask(:,:,1) * zconv )  
     189         &              ) * e1e2t ) * zconv 
    186190      ! salt flux 
    187       zsfx  = glob_sum( ( sfx + diag_smvi ) * e1e2t * tmask(:,:,1) * zconv ) * rday 
     191      zsfx  = glob_sum( ( sfx + diag_smvi ) * e1e2t ) * zconv * rday 
    188192      ! water flux 
    189       zvfx  = glob_sum( ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t * tmask(:,:,1) * zconv ) * rday 
     193      zvfx  = glob_sum( ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) * zconv * rday 
    190194 
    191195      ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice)  
    192       zarea   = glob_sum( SUM( a_i + epsi10, dim=3 ) * e1e2t * zconv ) ! in 1.e9 m2 
     196      zarea   = glob_sum( SUM( a_i + epsi10, dim=3 ) * e1e2t ) * zconv ! in 1.e9 m2 
    193197      zv_sill = zarea * 2.5e-5 
    194198      zs_sill = zarea * 25.e-5 
    195199      zh_sill = zarea * 10.e-5 
    196200 
    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 ) > zh_sill ) WRITE(numout,*) 'violation hfx    [GW]           (',cd_routine,')  = ',(zhfx) 
    200  
     201      IF(lwp) THEN 
     202         IF( ABS( zvfx ) > zv_sill )   WRITE(numout,*) 'violation vfx    [Mt/day]       (',cd_routine,')  = ',(zvfx) 
     203         IF( ABS( zsfx ) > zs_sill )   WRITE(numout,*) 'violation sfx    [psu*Mt/day]   (',cd_routine,')  = ',(zsfx) 
     204         IF( ABS( zhfx ) > zh_sill )   WRITE(numout,*) 'violation hfx    [GW]           (',cd_routine,')  = ',(zhfx) 
     205      ENDIF 
     206      ! 
    201207   END SUBROUTINE ice_cons_final 
    202208 
    … …  
    671677   !!   Default option         Empty Module               No LIM3 sea-ice model 
    672678   !!-------------------------------------------------------------------------- 
    673 CONTAINS 
    674    SUBROUTINE ice_ctl     ! Empty routine 
    675    END SUBROUTINE ice_ctl 
    676    SUBROUTINE ice_prt     ! Empty routine 
    677    END SUBROUTINE ice_prt 
    678    SUBROUTINE ice_prt3D   ! Empty routine 
    679    END SUBROUTINE ice_prt3D 
    680679#endif 
     680 
    681681   !!====================================================================== 
    682682END MODULE icectl 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedia.F90

    r8426 r8486  
    11MODULE icedia 
    22   !!====================================================================== 
    3    !!                       ***  MODULE limdia_hsb   *** 
    4    !!  LIM-3 sea ice model :   diagnostics of ice model  
     3   !!                       ***  MODULE icedia  *** 
     4   !!  Sea-Ice model :   global budgets  
    55   !!====================================================================== 
    66   !! History :  3.4  ! 2012-10  (C. Rousset)  original code 
     7   !!            4.0  ! 2017-08  (C. Rousset)  fits nemo4.0 standards 
    78   !!---------------------------------------------------------------------- 
    89#if defined key_lim3 
    … …  
    1011   !!   'key_lim3'                                       LIM3 sea-ice model 
    1112   !!---------------------------------------------------------------------- 
    12    !!   lim_dia_hsb        : computation and output of the time evolution of keys variables 
    13    !!   lim_dia_hsb_init   : initialization and namelist read 
    14    !!---------------------------------------------------------------------- 
    15    USE ice             ! LIM-3: sea-ice variable 
    16    USE dom_oce         ! ocean domain 
    17    USE sbc_oce, ONLY: sfx         ! surface boundary condition: ocean fields 
    18    USE daymod          ! model calendar 
    19    USE phycst          ! physical constant 
    20    USE in_out_manager  ! I/O manager 
    21    USE lib_mpp         ! MPP library 
    22    USE timing          ! preformance summary 
    23    USE iom             ! I/O manager 
    24    USE lib_fortran     ! glob_sum 
    25    USE icerst          ! ice restart 
     13   !!    ice_dia      : diagnostic of the sea-ice global heat content, salt content and volume conservation 
     14   !!    ice_dia_init : initialization of budget calculation 
     15   !!    ice_dia_rst  : read/write budgets restart 
     16   !!---------------------------------------------------------------------- 
     17   USE ice            ! LIM-3: sea-ice variable 
     18   USE dom_oce        ! ocean domain 
     19   USE phycst         ! physical constant 
     20   USE daymod         ! model calendar 
     21   USE sbc_oce , ONLY : sfx   ! surface boundary condition: ocean fields 
     22   USE icerst         ! ice restart 
     23   ! 
     24   USE in_out_manager ! I/O manager 
     25   USE lib_mpp        ! MPP library 
     26   USE timing         ! preformance summary 
     27   USE iom            ! I/O manager 
     28   USE lib_fortran    ! glob_sum 
    2629 
    2730   IMPLICIT NONE 
    … …  
    3639   !! * Substitutions 
    3740#  include "vectopt_loop_substitute.h90" 
    38  
    39    !!---------------------------------------------------------------------- 
    40    !! NEMO/OPA 3.4 , NEMO Consortium (2012) 
     41   !!---------------------------------------------------------------------- 
     42   !! NEMO/ICE 4.0 , NEMO Consortium (2017) 
    4143   !! $Id: icedia.F90 8413 2017-08-07 17:05:39Z clem $ 
    4244   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4345   !!---------------------------------------------------------------------- 
    44  
    4546CONTAINS 
    4647 
    … …  
    4950      !!                  ***  ROUTINE ice_dia  *** 
    5051      !!      
    51       !! ** Purpose: Compute the ice global heat content, salt content and volume conservation 
    52       !!  
    53       !!--------------------------------------------------------------------------- 
    54       INTEGER, INTENT(in) :: kt    ! number of iteration 
    55       !! 
    56       real(wp)   ::   zbg_ivol, zbg_svol, zbg_area, zbg_isal, zbg_item ,zbg_stem 
    57       REAL(wp)   ::   z_frc_voltop, z_frc_volbot, z_frc_sal, z_frc_temtop, z_frc_tembot   
     52      !! ** Purpose:   Compute the sea-ice global heat content, salt content  
     53      !!             and volume conservation 
     54      !!--------------------------------------------------------------------------- 
     55      INTEGER, INTENT(in) ::   kt   ! ocean time step  
     56      !! 
     57      REAL(wp)   ::   zbg_ivol, zbg_item, zbg_area, zbg_isal 
     58      REAL(wp)   ::   zbg_svol, zbg_stem 
     59      REAL(wp)   ::   z_frc_voltop, z_frc_temtop, z_frc_sal 
     60      REAL(wp)   ::   z_frc_volbot, z_frc_tembot   
    5861      REAL(wp)   ::   zdiff_vol, zdiff_sal, zdiff_tem   
    5962      !!--------------------------------------------------------------------------- 
    … …  
    6265      IF( kt == nit000 .AND. lwp ) THEN 
    6366         WRITE(numout,*) 
    64          WRITE(numout,*)'icedia' 
     67         WRITE(numout,*)'icedia : outpout ice diagnostics (integrated over the domain)' 
    6568         WRITE(numout,*)'~~~~~~' 
    6669      ENDIF 
    6770 
     71!!gm glob_sum includes a " * tmask_i ", so remove  " * tmask(:,:,1) " 
     72 
    6873      ! ----------------------- ! 
    6974      ! 1 -  Contents ! 
    7075      ! ----------------------- ! 
    71       zbg_ivol = glob_sum( vt_i(:,:) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 )                  ! ice volume (km3) 
    72       zbg_svol = glob_sum( vt_s(:,:) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 )                  ! snow volume (km3) 
    73       zbg_area = glob_sum( at_i(:,:) * e1e2t(:,:) * tmask(:,:,1) * 1.e-6 )                  ! area (km2) 
    74       zbg_isal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 ) ! salt content (pss*km3) 
    75       zbg_item = glob_sum( et_i * e1e2t(:,:) * tmask(:,:,1) * 1.e-20 )                      ! heat content (1.e20 J) 
    76       zbg_stem = glob_sum( et_s * e1e2t(:,:) * tmask(:,:,1) * 1.e-20 )                      ! heat content (1.e20 J) 
     76      zbg_ivol = glob_sum( vt_i(:,:) * e1e2t(:,:) ) * 1.e-9                  ! ice volume (km3) 
     77      zbg_svol = glob_sum( vt_s(:,:) * e1e2t(:,:) ) * 1.e-9                  ! snow volume (km3) 
     78      zbg_area = glob_sum( at_i(:,:) * e1e2t(:,:) ) * 1.e-6                  ! area (km2) 
     79      zbg_isal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e1e2t(:,:) ) * 1.e-9 ! salt content (pss*km3) 
     80      zbg_item = glob_sum( et_i * e1e2t(:,:) ) * 1.e-20                      ! heat content (1.e20 J) 
     81      zbg_stem = glob_sum( et_s * e1e2t(:,:) ) * 1.e-20                      ! heat content (1.e20 J) 
    7782       
    7883      ! ---------------------------! 
    7984      ! 2 - Trends due to forcing  ! 
    8085      ! ---------------------------! 
    81       z_frc_volbot = r1_rau0 * glob_sum( - ( wfx_ice(:,:) + wfx_snw(:,:) + wfx_err_sub(:,:) ) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 )  ! freshwater flux ice/snow-ocean  
    82       z_frc_voltop = r1_rau0 * glob_sum( - ( wfx_sub(:,:) + wfx_spr(:,:) ) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 )                     ! freshwater flux ice/snow-atm 
    83       z_frc_sal    = r1_rau0 * glob_sum( - sfx(:,:) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 )                                            ! salt fluxes ice/snow-ocean 
    84       z_frc_tembot =           glob_sum( hfx_out(:,:) * e1e2t(:,:) * tmask(:,:,1) * 1.e-20 )                                         ! heat on top of ocean (and below ice) 
    85       z_frc_temtop =           glob_sum( hfx_in (:,:) * e1e2t(:,:) * tmask(:,:,1) * 1.e-20 )                                         ! heat on top of ice-coean 
     86      z_frc_volbot = r1_rau0 * glob_sum( - ( wfx_ice(:,:) + wfx_snw(:,:) + wfx_err_sub(:,:) ) * e1e2t(:,:) ) * 1.e-9  ! freshwater flux ice/snow-ocean  
     87      z_frc_voltop = r1_rau0 * glob_sum( - ( wfx_sub(:,:) + wfx_spr(:,:) ) * e1e2t(:,:) ) * 1.e-9                     ! freshwater flux ice/snow-atm 
     88      z_frc_sal    = r1_rau0 * glob_sum(   - sfx(:,:) * e1e2t(:,:) ) * 1.e-9                                          ! salt fluxes ice/snow-ocean 
     89      z_frc_tembot =           glob_sum( hfx_out(:,:) * e1e2t(:,:) ) * 1.e-20                                         ! heat on top of ocean (and below ice) 
     90      z_frc_temtop =           glob_sum( hfx_in (:,:) * e1e2t(:,:) ) * 1.e-20                                         ! heat on top of ice-coean 
    8691      ! 
    8792      frc_voltop  = frc_voltop  + z_frc_voltop  * rdt_ice ! km3 
    … …  
    9499      ! 3 -  Content variations ! 
    95100      ! ----------------------- ! 
    96       zdiff_vol = r1_rau0 * glob_sum( ( rhoic * vt_i(:,:) + rhosn * vt_s(:,:) - vol_loc_ini(:,:)  &  ! freshwater trend (km3)  
    97          &                            ) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 )  
    98       zdiff_sal = r1_rau0 * glob_sum( ( rhoic * SUM( smv_i(:,:,:), dim=3 ) - sal_loc_ini(:,:)     &  ! salt content trend (km3*pss) 
    99          &                            ) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 ) 
    100       zdiff_tem =           glob_sum( ( et_i(:,:) + et_s(:,:) - tem_loc_ini(:,:)                  &  ! heat content trend (1.e20 J) 
    101       !  &                            + SUM( qevap_ice * a_i_b, dim=3 ) &     !! clem: I think this line should be commented (but needs a check) 
    102          &                            ) * e1e2t(:,:) * tmask(:,:,1) * 1.e-20 ) 
     101      zdiff_vol = r1_rau0 * glob_sum( ( rhoic*vt_i(:,:) + rhosn*vt_s(:,:) - vol_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! freshwater trend (km3)  
     102      zdiff_sal = r1_rau0 * glob_sum( ( rhoic* SUM( smv_i(:,:,:), dim=3 ) - sal_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! salt content trend (km3*pss) 
     103      zdiff_tem =           glob_sum( ( et_i(:,:) + et_s(:,:)             - tem_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-20 ! heat content trend (1.e20 J) 
     104      !                               + SUM( qevap_ice * a_i_b, dim=3 )       !! clem: I think this term should not be there (but needs a check) 
    103105 
    104106      ! ----------------------- ! 
    … …  
    112114      ! 5 - Diagnostics writing ! 
    113115      ! ----------------------- ! 
    114       ! 
    115       IF( iom_use('ibgvolume') )  CALL iom_put( 'ibgvolume' , zdiff_vol        )   ! ice/snow volume  drift            (km3 equivalent ocean water)          
    116       IF( iom_use('ibgsaltco') )  CALL iom_put( 'ibgsaltco' , zdiff_sal        )   ! ice salt content drift            (psu*km3 equivalent ocean water) 
    117       IF( iom_use('ibgheatco') )  CALL iom_put( 'ibgheatco' , zdiff_tem        )   ! ice/snow heat content drift       (1.e20 J) 
    118       IF( iom_use('ibgheatfx') )  CALL iom_put( 'ibgheatfx' , zdiff_tem /      &   ! ice/snow heat flux drift          (W/m2) 
    119          &                                                    glob_sum( e1e2t(:,:) * tmask(:,:,1) * 1.e-20 * kt*rdt ) ) 
    120  
    121       IF( iom_use('ibgfrcvoltop') )  CALL iom_put( 'ibgfrcvoltop' , frc_voltop )   ! vol  forcing ice/snw-atm          (km3 equivalent ocean water)  
    122       IF( iom_use('ibgfrcvolbot') )  CALL iom_put( 'ibgfrcvolbot' , frc_volbot )   ! vol  forcing ice/snw-ocean        (km3 equivalent ocean water)  
    123       IF( iom_use('ibgfrcsal') )     CALL iom_put( 'ibgfrcsal'    , frc_sal    )   ! sal - forcing                     (psu*km3 equivalent ocean water)    
    124       IF( iom_use('ibgfrctemtop') )  CALL iom_put( 'ibgfrctemtop' , frc_temtop )   ! heat on top of ice/snw/ocean      (1.e20 J)    
    125       IF( iom_use('ibgfrctembot') )  CALL iom_put( 'ibgfrctembot' , frc_tembot )   ! heat on top of ocean(below ice)   (1.e20 J)    
    126       IF( iom_use('ibgfrchfxtop') )  CALL iom_put( 'ibgfrchfxtop' , frc_temtop / & ! heat on top of ice/snw/ocean      (W/m2)  
    127          &                                                    glob_sum( e1e2t(:,:) * tmask(:,:,1) * 1.e-20 * kt*rdt ) ) 
    128       IF( iom_use('ibgfrchfxbot') )  CALL iom_put( 'ibgfrchfxbot' , frc_tembot / & ! heat on top of ocean(below ice)   (W/m2)  
    129          &                                                    glob_sum( e1e2t(:,:) * tmask(:,:,1) * 1.e-20 * kt*rdt ) ) 
    130  
    131       IF( iom_use('ibgvol_tot' ) )  CALL iom_put( 'ibgvol_tot'  , zbg_ivol     )   ! ice volume                        (km3) 
    132       IF( iom_use('sbgvol_tot' ) )  CALL iom_put( 'sbgvol_tot'  , zbg_svol     )   ! snow volume                       (km3) 
    133       IF( iom_use('ibgarea_tot') )  CALL iom_put( 'ibgarea_tot' , zbg_area     )   ! ice area                          (km2) 
    134       IF( iom_use('ibgsalt_tot') )  CALL iom_put( 'ibgsalt_tot' , zbg_isal     )   ! ice salinity content              (pss*km3) 
    135       IF( iom_use('ibgheat_tot') )  CALL iom_put( 'ibgheat_tot' , zbg_item     )   ! ice heat content                  (1.e20 J) 
    136       IF( iom_use('sbgheat_tot') )  CALL iom_put( 'sbgheat_tot' , zbg_stem     )   ! snow heat content                 (1.e20 J) 
     116!!gm I don't understand the division by the ocean surface (i.e. glob_sum( e1e2t(:,:) ) * 1.e-20 * kt*rdt ) 
     117!!   and its multiplication bu kt ! is it really what we want ? what is this quantity ? 
     118!!   IF it is really what we want, compute it at kt=nit000, not 3 time by time-step ! 
     119!!   kt*rdt  : you mean rdtice ? 
     120!!gm 
     121      ! 
     122      IF( iom_use('ibgvolume')    )   CALL iom_put( 'ibgvolume' , zdiff_vol     )   ! ice/snow volume  drift            (km3 equivalent ocean water)          
     123      IF( iom_use('ibgsaltco')    )   CALL iom_put( 'ibgsaltco' , zdiff_sal     )   ! ice salt content drift            (psu*km3 equivalent ocean water) 
     124      IF( iom_use('ibgheatco')    )   CALL iom_put( 'ibgheatco' , zdiff_tem     )   ! ice/snow heat content drift       (1.e20 J) 
     125      IF( iom_use('ibgheatfx')    )   CALL iom_put( 'ibgheatfx' ,               &   ! ice/snow heat flux drift          (W/m2) 
     126         &                                                     zdiff_tem /glob_sum( e1e2t(:,:) * 1.e-20 * kt*rdt ) ) 
     127 
     128      IF( iom_use('ibgfrcvoltop') )   CALL iom_put( 'ibgfrcvoltop' , frc_voltop )   ! vol  forcing ice/snw-atm          (km3 equivalent ocean water)  
     129      IF( iom_use('ibgfrcvolbot') )   CALL iom_put( 'ibgfrcvolbot' , frc_volbot )   ! vol  forcing ice/snw-ocean        (km3 equivalent ocean water)  
     130      IF( iom_use('ibgfrcsal')    )   CALL iom_put( 'ibgfrcsal'    , frc_sal    )   ! sal - forcing                     (psu*km3 equivalent ocean water)    
     131      IF( iom_use('ibgfrctemtop') )   CALL iom_put( 'ibgfrctemtop' , frc_temtop )   ! heat on top of ice/snw/ocean      (1.e20 J)    
     132      IF( iom_use('ibgfrctembot') )   CALL iom_put( 'ibgfrctembot' , frc_tembot )   ! heat on top of ocean(below ice)   (1.e20 J)    
     133      IF( iom_use('ibgfrchfxtop') )   CALL iom_put( 'ibgfrchfxtop' ,            &   ! heat on top of ice/snw/ocean      (W/m2)  
     134         &                                                          frc_temtop / glob_sum( e1e2t(:,:) ) * 1.e-20 * kt*rdt  ) 
     135      IF( iom_use('ibgfrchfxbot') )   CALL iom_put( 'ibgfrchfxbot' ,            &   ! heat on top of ocean(below ice)   (W/m2)  
     136         &                                                          frc_tembot / glob_sum( e1e2t(:,:) ) * 1.e-20 * kt*rdt  ) 
     137 
     138      IF( iom_use('ibgvol_tot' )  )   CALL iom_put( 'ibgvol_tot'  , zbg_ivol     )   ! ice volume                       (km3) 
     139      IF( iom_use('sbgvol_tot' )  )   CALL iom_put( 'sbgvol_tot'  , zbg_svol     )   ! snow volume                      (km3) 
     140      IF( iom_use('ibgarea_tot')  )   CALL iom_put( 'ibgarea_tot' , zbg_area     )   ! ice area                         (km2) 
     141      IF( iom_use('ibgsalt_tot')  )   CALL iom_put( 'ibgsalt_tot' , zbg_isal     )   ! ice salinity content             (pss*km3) 
     142      IF( iom_use('ibgheat_tot')  )   CALL iom_put( 'ibgheat_tot' , zbg_item     )   ! ice heat content                 (1.e20 J) 
     143      IF( iom_use('sbgheat_tot')  )   CALL iom_put( 'sbgheat_tot' , zbg_stem     )   ! snow heat content                (1.e20 J) 
    137144      ! 
    138145      IF( lrst_ice )   CALL ice_dia_rst( 'WRITE', kt_ice ) 
    … …  
    174181         RETURN 
    175182      ENDIF 
    176  
     183      ! 
    177184      CALL ice_dia_rst( 'READ' )  !* read or initialize all required files 
    178185      ! 
    179186   END SUBROUTINE ice_dia_init 
    180187 
     188 
    181189   SUBROUTINE ice_dia_rst( cdrw, kt ) 
    182      !!--------------------------------------------------------------------- 
    183      !!                   ***  ROUTINE limdia_rst  *** 
    184      !!                      
    185      !! ** Purpose :   Read or write DIA file in restart file 
    186      !! 
    187      !! ** Method  :   use of IOM library 
    188      !!---------------------------------------------------------------------- 
    189      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
    190      INTEGER         , INTENT(in), OPTIONAL ::   kt     ! ice time-step 
    191      REAL(wp)                     ::   ziter 
    192      INTEGER                      ::   iter 
    193      ! 
    194      !!---------------------------------------------------------------------- 
    195      ! 
    196      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise  
    197         IF( ln_rstart ) THEN                   !* Read the restart file 
    198            ! 
    199            CALL iom_get( numrir, 'kt_ice' , ziter ) 
    200            IF(lwp) WRITE(numout,*) 
    201            IF(lwp) WRITE(numout,*) ' ice_dia_rst read at time step = ', ziter 
    202            IF(lwp) WRITE(numout,*) '~~~~~~~' 
    203            CALL iom_get( numrir, 'frc_voltop' , frc_voltop  ) 
    204            CALL iom_get( numrir, 'frc_volbot' , frc_volbot  ) 
    205            CALL iom_get( numrir, 'frc_temtop' , frc_temtop  ) 
    206            CALL iom_get( numrir, 'frc_tembot' , frc_tembot  ) 
    207            CALL iom_get( numrir, 'frc_sal'    , frc_sal     ) 
    208            CALL iom_get( numrir, jpdom_autoglo, 'vol_loc_ini', vol_loc_ini ) 
    209            CALL iom_get( numrir, jpdom_autoglo, 'tem_loc_ini', tem_loc_ini ) 
    210            CALL iom_get( numrir, jpdom_autoglo, 'sal_loc_ini', sal_loc_ini ) 
    211         ELSE 
    212            IF(lwp) WRITE(numout,*) 
    213            IF(lwp) WRITE(numout,*) ' ice_dia at initial state ' 
    214            IF(lwp) WRITE(numout,*) '~~~~~~~' 
    215            ! set trends to 0 
    216            frc_voltop  = 0._wp                                           
    217            frc_volbot  = 0._wp                                           
    218            frc_temtop  = 0._wp                                                  
    219            frc_tembot  = 0._wp                                                  
    220            frc_sal     = 0._wp                                                  
    221            ! record initial ice volume, salt and temp 
    222            vol_loc_ini(:,:) = rhoic * vt_i(:,:) + rhosn * vt_s(:,:)  ! ice/snow volume (kg/m2) 
    223            tem_loc_ini(:,:) = et_i(:,:) + et_s(:,:)                  ! ice/snow heat content (J) 
    224            sal_loc_ini(:,:) = rhoic * SUM( smv_i(:,:,:), dim=3 )     ! ice salt content (pss*kg/m2) 
    225             
    226        ENDIF 
    227  
    228      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file 
    229         !                                   ! ------------------- 
    230         iter = kt + nn_fsbc - 1   ! ice restarts are written at kt == nitrst - nn_fsbc + 1 
    231  
    232         IF( iter == nitrst ) THEN 
    233            IF(lwp) WRITE(numout,*) 
    234            IF(lwp) WRITE(numout,*) ' ice_dia_rst write at time step = ', kt 
    235            IF(lwp) WRITE(numout,*) '~~~~~~~' 
    236         ENDIF 
    237  
    238         ! Write in numriw (if iter == nitrst) 
    239         ! ------------------  
    240         CALL iom_rstput( iter, nitrst, numriw, 'frc_voltop' , frc_voltop  ) 
    241         CALL iom_rstput( iter, nitrst, numriw, 'frc_volbot' , frc_volbot  ) 
    242         CALL iom_rstput( iter, nitrst, numriw, 'frc_temtop' , frc_temtop  ) 
    243         CALL iom_rstput( iter, nitrst, numriw, 'frc_tembot' , frc_tembot  ) 
    244         CALL iom_rstput( iter, nitrst, numriw, 'frc_sal'    , frc_sal     ) 
    245         CALL iom_rstput( iter, nitrst, numriw, 'vol_loc_ini', vol_loc_ini ) 
    246         CALL iom_rstput( iter, nitrst, numriw, 'tem_loc_ini', tem_loc_ini ) 
    247         CALL iom_rstput( iter, nitrst, numriw, 'sal_loc_ini', sal_loc_ini ) 
    248         ! 
    249      ENDIF 
    250      ! 
     190      !!--------------------------------------------------------------------- 
     191      !!                   ***  ROUTINE limdia_rst  *** 
     192      !!                      
     193      !! ** Purpose :   Read or write DIA file in restart file 
     194      !! 
     195      !! ** Method  :   use of IOM library 
     196      !!---------------------------------------------------------------------- 
     197      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
     198      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step 
     199      ! 
     200      INTEGER  ::   iter    ! local integer 
     201      REAL(wp) ::   ziter   ! local scalar 
     202      !!---------------------------------------------------------------------- 
     203      ! 
     204      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise  
     205         IF( ln_rstart ) THEN                   !* Read the restart file 
     206            ! 
     207            CALL iom_get( numrir, 'kt_ice' , ziter ) 
     208            IF(lwp) WRITE(numout,*) 
     209            IF(lwp) WRITE(numout,*) ' ice_dia_rst read at time step = ', ziter 
     210            IF(lwp) WRITE(numout,*) '~~~~~~~' 
     211            CALL iom_get( numrir, 'frc_voltop' , frc_voltop  ) 
     212            CALL iom_get( numrir, 'frc_volbot' , frc_volbot  ) 
     213            CALL iom_get( numrir, 'frc_temtop' , frc_temtop  ) 
     214            CALL iom_get( numrir, 'frc_tembot' , frc_tembot  ) 
     215            CALL iom_get( numrir, 'frc_sal'    , frc_sal     ) 
     216            CALL iom_get( numrir, jpdom_autoglo, 'vol_loc_ini', vol_loc_ini ) 
     217            CALL iom_get( numrir, jpdom_autoglo, 'tem_loc_ini', tem_loc_ini ) 
     218            CALL iom_get( numrir, jpdom_autoglo, 'sal_loc_ini', sal_loc_ini ) 
     219         ELSE 
     220            IF(lwp) WRITE(numout,*) 
     221            IF(lwp) WRITE(numout,*) ' ice_dia at initial state ' 
     222            IF(lwp) WRITE(numout,*) '~~~~~~~' 
     223            ! set trends to 0 
     224            frc_voltop  = 0._wp                                           
     225            frc_volbot  = 0._wp                                           
     226            frc_temtop  = 0._wp                                                  
     227            frc_tembot  = 0._wp                                                  
     228            frc_sal     = 0._wp                                                  
     229            ! record initial ice volume, salt and temp 
     230            vol_loc_ini(:,:) = rhoic * vt_i(:,:) + rhosn * vt_s(:,:)  ! ice/snow volume (kg/m2) 
     231            tem_loc_ini(:,:) = et_i(:,:) + et_s(:,:)                  ! ice/snow heat content (J) 
     232            sal_loc_ini(:,:) = rhoic * SUM( smv_i(:,:,:), dim=3 )     ! ice salt content (pss*kg/m2) 
     233         ENDIF 
     234         ! 
     235      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file 
     236         !                                   ! ------------------- 
     237         iter = kt + nn_fsbc - 1   ! ice restarts are written at kt == nitrst - nn_fsbc + 1 
     238         ! 
     239         IF( iter == nitrst ) THEN 
     240            IF(lwp) WRITE(numout,*) 
     241            IF(lwp) WRITE(numout,*) ' ice_dia_rst write at time step = ', kt 
     242            IF(lwp) WRITE(numout,*) '~~~~~~~' 
     243         ENDIF 
     244         ! 
     245         ! Write in numriw (if iter == nitrst) 
     246         ! ------------------  
     247         CALL iom_rstput( iter, nitrst, numriw, 'frc_voltop' , frc_voltop  ) 
     248         CALL iom_rstput( iter, nitrst, numriw, 'frc_volbot' , frc_volbot  ) 
     249         CALL iom_rstput( iter, nitrst, numriw, 'frc_temtop' , frc_temtop  ) 
     250         CALL iom_rstput( iter, nitrst, numriw, 'frc_tembot' , frc_tembot  ) 
     251         CALL iom_rstput( iter, nitrst, numriw, 'frc_sal'    , frc_sal     ) 
     252         CALL iom_rstput( iter, nitrst, numriw, 'vol_loc_ini', vol_loc_ini ) 
     253         CALL iom_rstput( iter, nitrst, numriw, 'tem_loc_ini', tem_loc_ini ) 
     254         CALL iom_rstput( iter, nitrst, numriw, 'sal_loc_ini', sal_loc_ini ) 
     255         ! 
     256      ENDIF 
     257      ! 
    251258   END SUBROUTINE ice_dia_rst 
    252259  
    … …  
    255262   !!   Default option :         Empty module          NO LIM sea-ice model 
    256263   !!---------------------------------------------------------------------- 
    257 CONTAINS 
    258    SUBROUTINE ice_dia          ! Empty routine 
    259    END SUBROUTINE ice_dia 
    260264#endif 
     265 
    261266   !!====================================================================== 
    262267END MODULE icedia 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90

    r8426 r8486  
    22   !!====================================================================== 
    33   !!                       ***  MODULE  iceforcing  *** 
    4    !! call surface forcing for sea ice model 
     4   !! Sea-Ice :   air-ice forcing fields 
    55   !!===================================================================== 
    66   !! History :  4.0  ! 2017-08  (C. Rousset) Original code 
    … …  
    1010   !!   'key_lim3' :                                  LIM 3.0 sea-ice model 
    1111   !!---------------------------------------------------------------------- 
    12    USE oce             ! ocean dynamics and tracers 
    13    USE dom_oce         ! ocean space and time domain 
    14    USE ice             ! LIM-3: ice variables 
     12   USE oce            ! ocean dynamics and tracers 
     13   USE dom_oce        ! ocean space and time domain 
     14   USE ice            ! sea-ice variables 
     15   USE sbc_oce        ! Surface boundary condition: ocean fields 
     16   USE sbc_ice        ! Surface boundary condition: ice   fields 
     17   USE usrdef_sbc     ! user defined: surface boundary condition 
     18   USE sbcblk         ! Surface boundary condition: bulk 
     19   USE sbccpl         ! Surface boundary condition: coupled interface 
     20   USE icealb         ! ice albedo 
    1521   ! 
    16    USE sbc_oce         ! Surface boundary condition: ocean fields 
    17    USE sbc_ice         ! Surface boundary condition: ice   fields 
    18    USE usrdef_sbc      ! user defined: surface boundary condition 
    19    USE sbcblk          ! Surface boundary condition: bulk 
    20    USE sbccpl          ! Surface boundary condition: coupled interface 
    21    USE icealb          ! ice albedo 
    22    ! 
    23    USE iom             ! I/O manager library 
    24    USE in_out_manager  ! I/O manager 
    25    USE lbclnk          ! lateral boundary condition - MPP link 
    26    USE lib_mpp         ! MPP library 
    27    USE lib_fortran     ! 
    28    USE timing          ! Timing 
     22   USE iom            ! I/O manager library 
     23   USE in_out_manager ! I/O manager 
     24   USE lbclnk         ! lateral boundary condition - MPP link 
     25   USE lib_mpp        ! MPP library 
     26   USE lib_fortran    ! 
     27   USE timing         ! Timing 
    2928 
    3029   IMPLICIT NONE 
    … …  
    3736#  include "vectopt_loop_substitute.h90" 
    3837   !!---------------------------------------------------------------------- 
    39    !! NEMO/OPA 4.0 , UCL NEMO Consortium (2011) 
     38   !! NEMO/ICE 4.0 , UCL NEMO Consortium (2017) 
    4039   !! $Id: icestp.F90 8319 2017-07-11 15:00:44Z clem $ 
    4140   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    113112      !!                alb_ice                                  = albedo above sea ice 
    114113      !!--------------------------------------------------------------------- 
    115       INTEGER, INTENT(in) ::   kt      ! ocean time step 
    116       INTEGER, INTENT(in) ::   ksbc    ! type of sbc flux ( 1 = user defined formulation,  
    117                                        !                    3 = bulk formulation, 
    118                                        !                    4 = Pure Coupled formulation) 
     114      INTEGER, INTENT(in) ::   kt     ! ocean time step 
     115      INTEGER, INTENT(in) ::   ksbc   ! flux formulation (user defined, bulk or Pure Coupled) 
     116      ! 
    119117      INTEGER  ::   ji, jj, jl                                ! dummy loop index 
    120118      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_os, zalb_cs  ! ice albedo under overcast/clear sky 
    121119      REAL(wp), DIMENSION(jpi,jpj)     ::   zalb              ! 2D workspace 
    122120      !!---------------------------------------------------------------------- 
    123  
     121      ! 
    124122      IF( nn_timing == 1 )   CALL timing_start('ice_forcing_flx') 
    125123 
    … …  
    136134      DO jl = 1, jpl 
    137135         DO jj = 2, jpjm1 
    138             DO ji = 2, jpim1                                 
     136            DO ji = 2, jpim1 
     137!!gm cldf_ice is a real, DOCTOR naming rule: start with cd means CHARACTER passed in argument ! 
    139138               alb_ice(ji,jj,jl) = ( 1. - cldf_ice ) * zalb_cs(ji,jj,jl) + cldf_ice * zalb_os(ji,jj,jl) 
    140139            END DO 
    … …  
    143142      CALL lbc_lnk( alb_ice, 'T', 1. ) 
    144143       
    145       ! --- fluxes over sea ice --- ! 
    146       SELECT CASE( ksbc ) 
    147          CASE( jp_usr )                                          ! user defined formulation 
    148                                    CALL usrdef_sbc_ice_flx( kt ) 
    149  
    150          CASE( jp_blk )                                          ! bulk formulation 
    151                                    CALL blk_ice_flx( t_su, alb_ice ) 
    152             IF( ln_mixcpl      )   CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su ) 
    153             IF( nn_limflx /= 2 )   CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
    154  
    155          CASE ( jp_purecpl )                                     ! coupled formulation 
    156                                    CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su ) 
    157             IF( nn_limflx == 2 )   CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
    158       END SELECT 
    159  
    160       ! --- albedo output --- ! 
    161       WHERE( at_i_b <= epsi06 )  ;  zalb(:,:) = rn_alb_oce 
    162       ELSEWHERE                  ;  zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) / at_i_b 
    163       END WHERE 
    164       IF( iom_use('icealb'  ) )  CALL iom_put( "icealb" , zalb(:,:) )   ! ice albedo 
    165  
    166       zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) + rn_alb_oce * ( 1._wp - at_i_b ) 
    167       IF( iom_use('albedo'  ) )  CALL iom_put( "albedo" , zalb(:,:) )   ! surface albedo 
    168       ! 
     144      SELECT CASE( ksbc )      !==  fluxes over sea ice  ==! 
     145      ! 
     146      CASE( jp_usr )                !--- user defined formulation 
     147                                CALL usrdef_sbc_ice_flx( kt ) 
     148         ! 
     149      CASE( jp_blk )                !--- bulk formulation 
     150                                CALL blk_ice_flx( t_su, alb_ice )    !  
     151         IF( ln_mixcpl      )   CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su ) 
     152         IF( nn_limflx /= 2 )   CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
     153         ! 
     154      CASE ( jp_purecpl )           !--- coupled formulation 
     155                                CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su ) 
     156         IF( nn_limflx == 2 )   CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
     157         ! 
     158      END SELECT 
     159 
     160      IF( iom_use('icealb') ) THEN    !--- output ice albedo 
     161         WHERE( at_i_b <= epsi06 )   ;   zalb(:,:) = rn_alb_oce 
     162         ELSEWHERE                   ;   zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) / at_i_b 
     163         END WHERE 
     164         CALL iom_put( "icealb" , zalb(:,:) )   ! ice albedo 
     165      ENDIF 
     166 
     167      IF( iom_use('albedo') ) THEN  !--- surface albedo 
     168         zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) + rn_alb_oce * ( 1._wp - at_i_b ) 
     169         CALL iom_put( "albedo" , zalb(:,:) ) 
     170      ENDIF 
    169171      ! 
    170172      IF( nn_timing == 1 )   CALL timing_stop('ice_forcing_flx') 
    171  
     173      ! 
    172174   END SUBROUTINE ice_forcing_flx 
    173175 
    … …  
    177179      !!                  ***  ROUTINE ice_flx_dist  *** 
    178180      !! 
    179       !! ** Purpose :   update the ice surface boundary condition by averaging and / or 
    180       !!                redistributing fluxes on ice categories 
     181      !! ** Purpose :   update the ice surface boundary condition by averaging  
     182      !!              and/or redistributing fluxes on ice categories 
    181183      !! 
    182184      !! ** Method  :   average then redistribute 
    … …  
    208210      IF( nn_timing == 1 )  CALL timing_start('ice_flx_dist') 
    209211      ! 
    210       SELECT CASE( k_limflx )                              !==  averaged on all ice categories  ==! 
     212      SELECT CASE( k_limflx )       !==  averaged on all ice categories  ==! 
     213      ! 
    211214      CASE( 0 , 1 ) 
    212          ! 
    213          z_qns_m  (:,:) = fice_ice_ave ( pqns_ice (:,:,:) ) 
    214          z_qsr_m  (:,:) = fice_ice_ave ( pqsr_ice (:,:,:) ) 
    215          z_dqn_m  (:,:) = fice_ice_ave ( pdqn_ice (:,:,:) ) 
    216          z_evap_m (:,:) = fice_ice_ave ( pevap_ice (:,:,:) ) 
    217          z_devap_m(:,:) = fice_ice_ave ( pdevap_ice (:,:,:) ) 
     215         z_qns_m  (:,:) = fice_ice_ave( pqns_ice  (:,:,:) ) 
     216         z_qsr_m  (:,:) = fice_ice_ave( pqsr_ice  (:,:,:) ) 
     217         z_dqn_m  (:,:) = fice_ice_ave( pdqn_ice  (:,:,:) ) 
     218         z_evap_m (:,:) = fice_ice_ave( pevap_ice (:,:,:) ) 
     219         z_devap_m(:,:) = fice_ice_ave( pdevap_ice(:,:,:) ) 
     220!!gm faster coding 
     221!    REAL(wp), DIMENSION(jpi,jpj) ::   z1_at_i   !  
     222! ... 
     223!      WHERE ( at_i (:,:) > 0._wp )   ; z1_at_i(:,:) = 1._wp / at_i (:,:) 
     224!      ELSEWHERE                      ; z1_at_i(:,:) = 0._wp 
     225!      END WHERE 
     226!      z_qns_m  (:,:) = SUM( a_i(:,:,:) * pqns_ice  (:,:,:) , dim=3 ) * z1_at_i(:,:) 
     227!      z_qsr_m  (:,:) = SUM( a_i(:,:,:) * pqsr_ice  (:,:,:) , dim=3 ) * z1_at_i(:,:) 
     228!      z_dqn_m  (:,:) = SUM( a_i(:,:,:) * pdqn_ice  (:,:,:) , dim=3 ) * z1_at_i(:,:) 
     229!      z_evap_m (:,:) = SUM( a_i(:,:,:) * pevap_ice (:,:,:) , dim=3 ) * z1_at_i(:,:) 
     230!      z_devap_m(:,:) = SUM( a_i(:,:,:) * pdevap_ice(:,:,:) , dim=3 ) * z1_at_i(:,:) 
     231!! and remove the 2 functions: fice_ice_ave and fice_cell_ave 
     232!!gm 
    218233         DO jl = 1, jpl 
    219             pdqn_ice  (:,:,jl) = z_dqn_m(:,:) 
     234            pdqn_ice  (:,:,jl) = z_dqn_m  (:,:) 
    220235            pdevap_ice(:,:,jl) = z_devap_m(:,:) 
    221          END DO 
    222          ! 
    223          DO jl = 1, jpl 
    224             pqns_ice (:,:,jl) = z_qns_m(:,:) 
    225             pqsr_ice (:,:,jl) = z_qsr_m(:,:) 
    226             pevap_ice(:,:,jl) = z_evap_m(:,:) 
    227          END DO 
    228          ! 
    229       END SELECT 
    230       ! 
    231       SELECT CASE( k_limflx )                              !==  redistribution on all ice categories  ==! 
     236            pqns_ice  (:,:,jl) = z_qns_m (:,:) 
     237            pqsr_ice  (:,:,jl) = z_qsr_m (:,:) 
     238            pevap_ice (:,:,jl) = z_evap_m(:,:) 
     239         END DO 
     240         ! 
     241      END SELECT 
     242      ! 
     243      SELECT CASE( k_limflx )       !==  redistribution on all ice categories  ==! 
    232244      CASE( 1 , 2 ) 
    233245         ! 
    234          zalb_m(:,:) = fice_ice_ave ( palb_ice (:,:,:) ) 
    235          ztem_m(:,:) = fice_ice_ave ( ptn_ice (:,:,:) ) 
     246         zalb_m(:,:) = fice_ice_ave( palb_ice(:,:,:) ) 
     247         ztem_m(:,:) = fice_ice_ave( ptn_ice (:,:,:) ) 
    236248         DO jl = 1, jpl 
    237249            pqns_ice (:,:,jl) = pqns_ice (:,:,jl) + pdqn_ice  (:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) ) 
    … …  
    246258   END SUBROUTINE ice_flx_dist 
    247259 
    248  
     260!!gm TO BE REMOVED ====>>>>> 
    249261   FUNCTION fice_cell_ave ( ptab ) 
    250262      !!-------------------------------------------------------------------------- 
    251263      !! * Compute average over categories, for grid cell (ice covered and free ocean) 
    252264      !!-------------------------------------------------------------------------- 
    253       REAL(wp), DIMENSION(jpi,jpj) :: fice_cell_ave 
    254       REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT (in) :: ptab 
     265      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(in) :: ptab 
     266      REAL(wp), DIMENSION(jpi,jpj)                 :: fice_cell_ave 
    255267      INTEGER :: jl 
    256268      !!-------------------------------------------------------------------------- 
    257  
    258       fice_cell_ave (:,:) = 0._wp 
    259       DO jl = 1, jpl 
    260          fice_cell_ave (:,:) = fice_cell_ave (:,:) + a_i (:,:,jl) * ptab (:,:,jl) 
     269      fice_cell_ave(:,:) = a_i(:,:,1) * ptab (:,:,1) 
     270      DO jl = 2, jpl 
     271         fice_cell_ave(:,:) = fice_cell_ave(:,:) + a_i(:,:,jl) * ptab (:,:,jl) 
    261272      END DO 
    262  
    263273   END FUNCTION fice_cell_ave 
    264274 
    … …  
    268278      !! * Compute average over categories, for ice covered part of grid cell 
    269279      !!-------------------------------------------------------------------------- 
    270       REAL(wp), DIMENSION(jpi,jpj) :: fice_ice_ave 
    271       REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(in) :: ptab 
    272      !!-------------------------------------------------------------------------- 
    273  
     280      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(in) ::   ptab   ! 
     281      REAL(wp), DIMENSION(jpi,jpj)                 :: fice_ice_ave 
     282      !!-------------------------------------------------------------------------- 
    274283      WHERE ( at_i (:,:) > 0.0_wp ) ; fice_ice_ave (:,:) = fice_cell_ave( ptab (:,:,:) ) / at_i (:,:) 
    275284      ELSEWHERE                     ; fice_ice_ave (:,:) = 0.0_wp 
    276285      END WHERE 
    277        
    278286   END FUNCTION fice_ice_ave 
    279287 
     288!!gm <<<<<<====  end of TO BE REMOVED  
     289 
     290#else 
     291   !!---------------------------------------------------------------------- 
     292   !!   Default option :         Empty module          NO LIM sea-ice model 
     293   !!---------------------------------------------------------------------- 
    280294#endif 
    281295 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceist.F90

    r8426 r8486  
    55   !!====================================================================== 
    66   !! History :  2.0  ! 2004-01 (C. Ethe, G. Madec)  Original code 
    7    !!            3.0  ! 2011-02 (G. Madec) dynamical allocation 
    8    !!             -   ! 2014    (C. Rousset) add N/S initializations 
     7   !!            3.0  ! 2007    (M. Vancoppenolle)  Rewrite for ice cats 
     8   !!            3.0  ! 2009-11 (M. Vancoppenolle)  Enhanced version for ice cats 
     9   !!            3.0  ! 2011-02 (G. Madec)  dynamical allocation 
     10   !!             -   ! 2014    (C. Rousset)  add N/S initializations 
    911   !!---------------------------------------------------------------------- 
    1012#if defined key_lim3 
    1113   !!---------------------------------------------------------------------- 
    12    !!   'key_lim3' :                                    LIM3 sea-ice model 
    13    !!---------------------------------------------------------------------- 
    14    !!   ice_ist      :  Initialisation of diagnostics ice variables 
    15    !!   ice_ist_init :  initialization of ice state and namelist read 
     14   !!   'key_lim3'                                       LIM3 sea-ice model 
     15   !!---------------------------------------------------------------------- 
     16   !!   ice_ist       :  initialization of diagnostics ice variables 
     17   !!   ice_ist_init  :  initialization of ice state and namelist read 
    1618   !!---------------------------------------------------------------------- 
    1719   USE phycst         ! physical constant 
    … …  
    3436   PRIVATE 
    3537 
    36    PUBLIC   ice_ist      ! routine called by lim_init.F90 
     38   PUBLIC   ice_ist      ! called by icestp.F90 
    3739 
    3840   INTEGER , PARAMETER ::   jpfldi = 6           ! maximum number of files to read 
    … …  
    4547   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si  ! structure of input fields (file informations, fields read) 
    4648   !!---------------------------------------------------------------------- 
    47    !!   LIM 3.0,  UCL-LOCEAN-IPSL (2008) 
     49   !! NEMO/ICE 4.0 , NEMO Consortium (2017) 
    4850   !! $Id: iceist.F90 8378 2017-07-26 13:55:59Z clem $ 
    4951   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    … …  
    5153CONTAINS 
    5254 
     55!!gm  better name:  ice_istate 
    5356   SUBROUTINE ice_ist 
    5457      !!------------------------------------------------------------------- 
    … …  
    7174      !! ** Notes   : o_i, t_su, t_s, t_i, s_i must be filled everywhere, even 
    7275      !!              where there is no ice (clem: I do not know why, is it mandatory?)  
    73       !! 
    74       !! History : 
    75       !!   2.0  !  01-04  (C. Ethe, G. Madec)  Original code 
    76       !!   3.0  !  2007   (M. Vancoppenolle)   Rewrite for ice cats 
    77       !!   4.0  !  09-11  (M. Vancoppenolle)   Enhanced version for ice cats 
    7876      !!-------------------------------------------------------------------- 
    79       !! * Local variables 
    8077      INTEGER    :: ji, jj, jk, jl             ! dummy loop indices 
    8178      REAL(wp)   :: ztmelts, zdh 
    … …  
    591588   !!   Default option :         Empty module          NO LIM sea-ice model 
    592589   !!---------------------------------------------------------------------- 
    593 CONTAINS 
    594    SUBROUTINE ice_ist          ! Empty routine 
    595    END SUBROUTINE ice_ist 
    596590#endif 
    597591 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90

    r8426 r8486  
    1111#if defined key_lim3 
    1212   !!---------------------------------------------------------------------- 
    13    !!   'key_lim3' :                                   LIM3 sea-ice model 
     13   !!   'key_lim3'                                       LIM3 sea-ice model 
    1414   !!---------------------------------------------------------------------- 
    1515   !!   ice_itd_rem   : 
    … …  
    1818   !!   ice_itd_shiftice : 
    1919   !!---------------------------------------------------------------------- 
    20    USE par_oce          ! ocean parameters 
    21    USE dom_oce          ! ocean domain 
    22    USE phycst           ! physical constants (ocean directory)  
    23    USE ice1D            ! LIM-3 thermodynamic variables 
    24    USE ice              ! LIM-3 variables 
    25    USE icectl           ! conservation tests 
    26    USE icetab 
     20   USE par_oce        ! ocean parameters 
     21   USE dom_oce        ! ocean domain 
     22   USE phycst         ! physical constants  
     23   USE ice1D          ! sea-ice: thermodynamic variables 
     24   USE ice            ! sea-ice: variables 
     25   USE icectl         ! sea-ice: conservation tests 
     26   USE icetab         ! sea-ice: convert 1D<=>2D 
    2727   ! 
    28    USE prtctl           ! Print control 
    29    USE in_out_manager   ! I/O manager 
    30    USE lib_mpp          ! MPP library 
    31    USE lib_fortran      ! to use key_nosignedzero 
     28   USE prtctl         ! Print control 
     29   USE in_out_manager ! I/O manager 
     30   USE lib_mpp        ! MPP library 
     31   USE lib_fortran    ! to use key_nosignedzero 
    3232 
    3333   IMPLICIT NONE 
    … …  
    3838 
    3939   !!---------------------------------------------------------------------- 
    40    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 
     40   !! NEMO/ICE 4.0 , NEMO Consortium (2017) 
    4141   !! $Id: iceitd.F90 8420 2017-08-08 12:18:46Z clem $ 
    4242   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    8686      !  1) Identify grid cells with ice 
    8787      !----------------------------------------------------------------------------------------------- 
    88       nidx = 0 ; idxice(:) = 0 
     88      nidx = 0   ;  idxice(:) = 0 
    8989      DO jj = 1, jpj 
    9090         DO ji = 1, jpi 
    … …  
    118118         ! --- New boundaries for category 1:jpl-1 --- ! 
    119119         DO jl = 1, jpl - 1 
    120  
     120            ! 
    121121            DO ji = 1, nidx 
    122122               ! 
    … …  
    134134                  zhbnew(ji,jl) = hi_max(jl) 
    135135               ENDIF 
    136              
     136               ! 
    137137               ! --- 2 conditions for remapping --- ! 
    138138               ! 1) hn(t+1)+espi < Hn* < hn+1(t+1)-epsi                
    … …  
    148148            END DO 
    149149         END DO 
    150           
     150         ! 
    151151         ! --- New boundaries for category jpl --- ! 
    152152         DO ji = 1, nidx 
    … …  
    164164            IF( ht_ib_2d(ji,1) > ( hi_max(1) - epsi10 ) )   idxice(ji) = 0 
    165165         END DO 
    166  
     166         ! 
    167167         !----------------------------------------------------------------------------------------------- 
    168168         !  3) Identify cells where remapping 
    … …  
    178178         idxice(:) = idxice2(:) 
    179179         nidx      = nidx2 
    180           
     180         ! 
    181181      ENDIF 
    182182    
    … …  
    186186      !----------------------------------------------------------------------------------------------- 
    187187      IF( nidx > 0 ) THEN 
    188  
    189          zhb0(:) = hi_max(0) ; zhb1(:) = hi_max(1) 
    190          g0(:,:) = 0._wp     ; g1(:,:) = 0._wp  
    191          hL(:,:) = 0._wp     ; hR(:,:) = 0._wp  
    192           
     188         ! 
     189         zhb0(:) = hi_max(0)   ;  zhb1(:) = hi_max(1) 
     190         g0(:,:) = 0._wp       ;  g1(:,:) = 0._wp  
     191         hL(:,:) = 0._wp       ;  hR(:,:) = 0._wp  
     192         ! 
    193193         DO jl = 1, jpl 
    194  
     194            ! 
    195195            CALL tab_2d_1d( nidx, idxice(1:nidx), ht_ib_1d(1:nidx), ht_i_b(:,:,jl) ) 
    196196            CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,jl)   ) 
    197197            CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl)    ) 
    198198            CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl)    ) 
    199  
     199            ! 
    200200            IF( jl == 1 ) THEN 
    201                 
     201               !   
    202202               ! --- g(h) for category 1 --- ! 
    203                CALL ice_itd_glinear( zhb0(1:nidx), zhb1(1:nidx), ht_ib_1d(1:nidx), a_i_1d(1:nidx),  &   ! in 
    204                   &                  g0(1:nidx,1), g1(1:nidx,1), hL(1:nidx,1)    , hR(1:nidx,1) )       ! out 
    205                 
     203               CALL ice_itd_glinear( zhb0(1:nidx)  , zhb1(1:nidx)  , ht_ib_1d(1:nidx)  , a_i_1d(1:nidx)  ,  &   ! in 
     204                  &                  g0  (1:nidx,1), g1  (1:nidx,1), hL      (1:nidx,1), hR    (1:nidx,1)   )   ! out 
     205                  ! 
    206206               ! Area lost due to melting of thin ice 
    207207               DO ji = 1, nidx 
    208                    
     208                  ! 
    209209                  IF( a_i_1d(ji) > epsi10 ) THEN 
    210                       
     210                     ! 
    211211                     zdh0 =  ht_i_1d(ji) - ht_ib_1d(ji)                 
    212212                     IF( zdh0 < 0.0 ) THEN      !remove area from category 1 
    … …  
    214214                        !Integrate g(1) from 0 to dh0 to estimate area melted 
    215215                        zetamax = MIN( zdh0, hR(ji,1) ) - hL(ji,1) 
    216                          
     216                        ! 
    217217                        IF( zetamax > 0.0 ) THEN 
    218218                           zx1    = zetamax 
    … …  
    227227                           v_i_1d(ji)  = a_i_1d(ji) * ht_i_1d(ji) ! clem-useless ? 
    228228                        ENDIF 
    229                          
     229                        ! 
    230230                     ELSE ! if ice accretion zdh0 > 0 
    231231                        ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1) 
    232232                        zhbnew(ji,0) = MIN( zdh0, hi_max(1) )  
    233233                     ENDIF 
    234                       
     234                     ! 
    235235                  ENDIF 
    236                    
     236                  ! 
    237237               END DO 
    238  
     238               ! 
    239239               CALL tab_1d_2d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,jl)   ) 
    240240               CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl)    ) 
    241241               CALL tab_1d_2d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl)    ) 
    242                 
     242               ! 
    243243            ENDIF ! jl=1 
    244              
     244            ! 
    245245            ! --- g(h) for each thickness category --- !   
    246             CALL ice_itd_glinear( zhbnew(1:nidx,jl-1), zhbnew(1:nidx,jl), ht_i_1d(1:nidx), a_i_1d(1:nidx),  &  ! in 
    247                &                  g0(1:nidx,jl)      , g1(1:nidx,jl)    , hL(1:nidx,jl)  , hR(1:nidx,jl) )     ! out 
    248  
     246            CALL ice_itd_glinear( zhbnew(1:nidx,jl-1), zhbnew(1:nidx,jl), ht_i_1d(1:nidx)   , a_i_1d(1:nidx)   ,  &   ! in 
     247               &                  g0    (1:nidx,jl  ), g1    (1:nidx,jl), hL     (1:nidx,jl), hR    (1:nidx,jl)   )   ! out 
     248            ! 
    249249         END DO 
    250250          
    … …  
    253253         !----------------------------------------------------------------------------------------------- 
    254254         DO jl = 1, jpl - 1 
    255              
     255            ! 
    256256            DO ji = 1, nidx 
    257                 
     257               ! 
    258258               ! left and right integration limits in eta space 
    259259               IF (zhbnew(ji,jl) > hi_max(jl)) THEN ! Hn* > Hn => transfer from jl to jl+1 
    … …  
    267267               ENDIF 
    268268               zetamax = MAX( zetamax, zetamin ) ! no transfer if etamax < etamin 
    269                 
     269               ! 
    270270               zx1  = zetamax - zetamin 
    271271               zwk1 = zetamin * zetamin 
    … …  
    278278               zdaice(ji,jl) = g1(ji,jcat)*zx2 + g0(ji,jcat)*zx1 
    279279               zdvice(ji,jl) = g1(ji,jcat)*zx3 + g0(ji,jcat)*zx2 + zdaice(ji,jl)*hL(ji,jcat) 
    280                 
     280               ! 
    281281            END DO 
    282282         END DO 
    … …  
    305305            ENDIF 
    306306         END DO 
    307  
     307         ! 
    308308         CALL tab_1d_2d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,1)   ) 
    309309         CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,1)    ) 
    310310         CALL tab_1d_2d( nidx, idxice(1:nidx), a_ip_1d (1:nidx), a_ip(:,:,1)    ) 
    311          
     311         ! 
    312312      ENDIF 
    313  
    314       IF( ln_limdiachk ) CALL ice_cons_hsm(1, 'iceitd_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    315  
     313      ! 
     314      IF( ln_limdiachk )   CALL ice_cons_hsm(1, 'iceitd_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     315      ! 
    316316   END SUBROUTINE ice_itd_rem 
    317317 
    … …  
    325325      !! ** Method  :   g(h) is linear and written as: g(eta) = g1(eta) + g0 
    326326      !!                with eta = h - HL 
    327       !!                 
    328327      !!------------------------------------------------------------------ 
    329328      REAL(wp), DIMENSION(:), INTENT(in   ) ::   HbL, HbR      ! left and right category boundaries 
    … …  
    333332      ! 
    334333      INTEGER  ::   ji           ! horizontal indices 
     334      REAL(wp) ::   z1_3 , z2_3  ! 1/3 , 2/3 
    335335      REAL(wp) ::   zh13         ! HbL + 1/3 * (HbR - HbL) 
    336336      REAL(wp) ::   zh23         ! HbL + 2/3 * (HbR - HbL) 
    … …  
    339339      !!------------------------------------------------------------------ 
    340340      ! 
    341          DO ji = 1, nidx 
    342             ! 
    343             IF( paice(ji) > epsi10  .AND. phice(ji) > 0._wp )  THEN 
    344  
    345                ! Initialize hL and hR 
    346                phL(ji) = HbL(ji) 
    347                phR(ji) = HbR(ji) 
    348  
    349                ! Change hL or hR if hice falls outside central third of range, 
    350                ! so that hice is in the central third of the range [HL HR] 
    351                zh13 = 1.0 / 3.0 * ( 2.0 * phL(ji) +       phR(ji) ) 
    352                zh23 = 1.0 / 3.0 * (       phL(ji) + 2.0 * phR(ji) ) 
    353  
    354                IF    ( phice(ji) < zh13 ) THEN   ;   phR(ji) = 3._wp * phice(ji) - 2._wp * phL(ji) ! move HR to the left 
    355                ELSEIF( phice(ji) > zh23 ) THEN   ;   phL(ji) = 3._wp * phice(ji) - 2._wp * phR(ji) ! move HL to the right 
    356                ENDIF 
    357  
    358                ! Compute coefficients of g(eta) = g0 + g1*eta 
    359                zdhr = 1._wp / (phR(ji) - phL(ji)) 
    360                zwk1 = 6._wp * paice(ji) * zdhr 
    361                zwk2 = ( phice(ji) - phL(ji) ) * zdhr 
    362                pg0(ji) = zwk1 * ( 2._wp / 3._wp - zwk2 )      ! Eq. 14 
    363                pg1(ji) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 ) ! Eq. 14 
    364                ! 
    365             ELSE  ! remap_flag = .false. or a_i < epsi10  
    366                phL(ji) = 0._wp 
    367                phR(ji) = 0._wp 
    368                pg0(ji) = 0._wp 
    369                pg1(ji) = 0._wp 
    370             ENDIF 
    371             ! 
    372          END DO 
     341      z1_3 = 1._wp / 3._wp 
     342      z2_3 = 2._wp / 3._wp 
     343      ! 
     344      DO ji = 1, nidx 
     345         ! 
     346         IF( paice(ji) > epsi10  .AND. phice(ji) > 0._wp )  THEN 
     347            ! 
     348            ! Initialize hL and hR 
     349            phL(ji) = HbL(ji) 
     350            phR(ji) = HbR(ji) 
     351            ! 
     352            ! Change hL or hR if hice falls outside central third of range, 
     353            ! so that hice is in the central third of the range [HL HR] 
     354            zh13 = z1_3 * ( 2._wp * phL(ji) +         phR(ji) ) 
     355            zh23 = z1_3 * (         phL(ji) + 2._wp * phR(ji) ) 
     356            ! 
     357            IF    ( phice(ji) < zh13 ) THEN   ;   phR(ji) = 3._wp * phice(ji) - 2._wp * phL(ji) ! move HR to the left 
     358            ELSEIF( phice(ji) > zh23 ) THEN   ;   phL(ji) = 3._wp * phice(ji) - 2._wp * phR(ji) ! move HL to the right 
     359            ENDIF 
     360            ! 
     361            ! Compute coefficients of g(eta) = g0 + g1*eta 
     362            zdhr = 1._wp / (phR(ji) - phL(ji)) 
     363            zwk1 = 6._wp * paice(ji) * zdhr 
     364            zwk2 = ( phice(ji) - phL(ji) ) * zdhr 
     365            pg0(ji) = zwk1 * ( z2_3 - zwk2 )                    ! Eq. 14 
     366            pg1(ji) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5_wp )   ! Eq. 14 
     367            ! 
     368         ELSE  ! remap_flag = .false. or a_i < epsi10  
     369            phL(ji) = 0._wp 
     370            phR(ji) = 0._wp 
     371            pg0(ji) = 0._wp 
     372            pg1(ji) = 0._wp 
     373         ENDIF 
     374         ! 
     375      END DO 
    373376      ! 
    374377   END SUBROUTINE ice_itd_glinear 
    … …  
    381384      !! ** Purpose :   shift ice across category boundaries, conserving everything 
    382385      !!              ( area, volume, energy, age*vol, and mass of salt ) 
    383       !! 
    384       !! ** Method  : 
    385386      !!------------------------------------------------------------------ 
    386387      INTEGER , DIMENSION(:,:), INTENT(in) ::   kdonor   ! donor category index 
    387388      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pdaice   ! ice area transferred across boundary 
    388389      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pdvice   ! ice volume transferred across boundary 
    389  
    390       INTEGER ::   ii,ij, ji, jj, jl, jl2, jl1, jk   ! dummy loop indices 
    391  
    392       REAL(wp), DIMENSION(jpij,jpl) ::   zaTsfn 
    393       REAL(wp), DIMENSION(jpij)     ::   zworka            ! temporary array used here 
    394       REAL(wp), DIMENSION(jpij)     ::   zworkv            ! temporary array used here 
    395  
     390      ! 
     391      INTEGER  ::   ji, jj, jl, jk     ! dummy loop indices 
     392      INTEGER  ::   ii, ij, jl2, jl1   ! local integers 
    396393      REAL(wp) ::   ztrans             ! ice/snow transferred 
     394      REAL(wp), DIMENSION(jpij)     ::   zworka, zworkv   ! workspace 
     395      REAL(wp), DIMENSION(jpij,jpl) ::   zaTsfn           !  -    - 
    397396      !!------------------------------------------------------------------ 
    398397          
    … …  
    421420      DO jl = 1, jpl - 1 
    422421         DO ji = 1, nidx 
    423              
     422            ! 
    424423            jl1 = kdonor(ji,jl) 
    425             IF    ( jl1 == jl  ) THEN   ;   jl2 = jl1+1 
    426             ELSE                        ;   jl2 = jl  
    427             ENDIF 
    428  
    429             rswitch    = MAX( 0._wp , SIGN( 1._wp , v_i_2d(ji,jl1) - epsi10 ) ) 
    430             zworkv(ji) = pdvice(ji,jl) / MAX( v_i_2d(ji,jl1), epsi10 ) * rswitch 
    431  
    432             rswitch    = MAX( 0._wp , SIGN( 1._wp , a_i_2d(ji,jl1) - epsi10 ) ) 
    433             zworka(ji) = pdaice(ji,jl) / MAX( a_i_2d(ji,jl1), epsi10 ) * rswitch         
    434                 
    435             ! Ice areas 
    436             a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl) 
     424            IF ( jl1 == jl  ) THEN   ;   jl2 = jl1+1 
     425            ELSE                     ;   jl2 = jl  
     426            ENDIF 
     427            ! 
     428!!gm use of rswitch not faster as there is already IF in the DO-loop ==>>> use IF which is more readable 
     429!            rswitch    = MAX( 0._wp , SIGN( 1._wp , v_i_2d(ji,jl1) - epsi10 ) ) 
     430!            zworkv(ji) = pdvice(ji,jl) / MAX( v_i_2d(ji,jl1), epsi10 ) * rswitch 
     431!            ! 
     432!            rswitch    = MAX( 0._wp , SIGN( 1._wp , a_i_2d(ji,jl1) - epsi10 ) ) 
     433!            zworka(ji) = pdaice(ji,jl) / MAX( a_i_2d(ji,jl1), epsi10 ) * rswitch         
     434 
     435            IF( v_i_2d(ji,jl1) >= epsi10 ) THEN   ;   zworkv(ji) = pdvice(ji,jl) / v_i_2d(ji,jl1) 
     436            ELSE                                  ;   zworkv(ji) = 0._wp 
     437            ENDIF 
     438            IF( a_i_2d(ji,jl1) >= epsi10 ) THEN   ;   zworka(ji) = pdaice(ji,jl) / a_i_2d(ji,jl1) 
     439            ELSE                                  ;   zworka(ji) = 0._wp 
     440            ENDIF 
     441!!gm end 
     442            ! 
     443            a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl)       ! Ice areas 
    437444            a_i_2d(ji,jl2) = a_i_2d(ji,jl2) + pdaice(ji,jl) 
    438                 
    439             ! Ice volumes 
    440             v_i_2d(ji,jl1) = v_i_2d(ji,jl1) - pdvice(ji,jl)  
     445            ! 
     446            v_i_2d(ji,jl1) = v_i_2d(ji,jl1) - pdvice(ji,jl)       ! Ice volumes 
    441447            v_i_2d(ji,jl2) = v_i_2d(ji,jl2) + pdvice(ji,jl) 
    442              
    443             ! Snow volumes 
    444             ztrans         = v_s_2d(ji,jl1) * zworkv(ji) 
     448            ! 
     449            ztrans         = v_s_2d(ji,jl1) * zworkv(ji)          ! Snow volumes 
    445450            v_s_2d(ji,jl1) = v_s_2d(ji,jl1) - ztrans 
    446451            v_s_2d(ji,jl2) = v_s_2d(ji,jl2) + ztrans  
    447                         
    448             ! Ice age  
    449             ztrans             = oa_i_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work 
     452            !           
     453            !                                                     ! Ice age  
     454            ztrans             = oa_i_2d(ji,jl1) * pdaice(ji,jl)  !!clem: should be * zworka(ji) but it does not work ???? 
    450455            oa_i_2d(ji,jl1)    = oa_i_2d(ji,jl1) - ztrans 
    451456            oa_i_2d(ji,jl2)    = oa_i_2d(ji,jl2) + ztrans 
    452              
    453             ! Ice salinity 
    454             ztrans             = smv_i_2d(ji,jl1) * zworkv(ji) 
     457            ! 
     458            ztrans             = smv_i_2d(ji,jl1) * zworkv(ji)    ! Ice salinity 
     459 
    455460            smv_i_2d(ji,jl1)   = smv_i_2d(ji,jl1) - ztrans 
    456461            smv_i_2d(ji,jl2)   = smv_i_2d(ji,jl2) + ztrans 
    457              
    458             ! Surface temperature 
    459             ztrans          = t_su_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work 
     462            ! 
     463            !                                                     ! Surface temperature 
     464            ztrans          = t_su_2d(ji,jl1) * pdaice(ji,jl)     !!clem: should be * zworka(ji) but it does not work ???? 
    460465            zaTsfn(ji,jl1)  = zaTsfn(ji,jl1) - ztrans 
    461466            zaTsfn(ji,jl2)  = zaTsfn(ji,jl2) + ztrans 
    462                 
     467            !   
    463468            ! MV MP 2016  
    464469            IF ( nn_pnd_scheme > 0 ) THEN 
    465                ! Pond fraction 
     470               !                                                  ! Pond fraction 
    466471               ztrans          = a_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work 
    467472               a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans 
    468473               a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans 
    469                 
    470                ! Pond volume (also proportional to da/a) 
     474               !                                                  ! Pond volume (also proportional to da/a) 
    471475               ztrans          = v_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work 
    472476               v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans 
    … …  
    474478            ENDIF 
    475479            ! END MV MP 2016  
    476                 
    477          END DO 
    478  
    479          ! Snow heat content 
    480          DO jk = 1, nlay_s 
    481  
     480         END DO 
     481         ! 
     482         DO jk = 1, nlay_s         !--- Snow heat content 
     483            ! 
    482484            DO ji = 1, nidx 
    483485               ii = MOD( idxice(ji) - 1, jpi ) + 1 
    484486               ij = ( idxice(ji) - 1 ) / jpi + 1 
    485  
     487               ! 
    486488               jl1 = kdonor(ji,jl) 
    487489               IF(jl1 == jl) THEN  ;  jl2 = jl+1 
    488490               ELSE                ;  jl2 = jl 
    489491               ENDIF 
    490                 
     492               ! 
    491493               ztrans            = e_s(ii,ij,jk,jl1) * zworkv(ji) 
    492494               e_s(ii,ij,jk,jl1) = e_s(ii,ij,jk,jl1) - ztrans 
    … …  
    495497         END DO 
    496498 
    497           
    498          ! Ice heat content 
    499          DO jk = 1, nlay_i 
     499         DO jk = 1, nlay_i         !--- Ice heat content 
    500500            DO ji = 1, nidx 
    501501               ii = MOD( idxice(ji) - 1, jpi ) + 1 
    502502               ij = ( idxice(ji) - 1 ) / jpi + 1 
    503                 
     503               ! 
    504504               jl1 = kdonor(ji,jl) 
    505505               IF(jl1 == jl) THEN  ;  jl2 = jl+1 
    506506               ELSE                ;  jl2 = jl 
    507507               ENDIF 
    508                 
     508               ! 
    509509               ztrans            = e_i(ii,ij,jk,jl1) * zworkv(ji) 
    510510               e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - ztrans 
    … …  
    512512            END DO 
    513513         END DO 
    514           
     514         ! 
    515515      END DO                   ! boundaries, 1 to jpl-1 
    516516       
    … …  
    521521         DO ji = 1, nidx 
    522522            IF ( a_i_2d(ji,jl) > epsi10 ) THEN  
    523                ht_i_2d(ji,jl)  =  v_i_2d   (ji,jl) / a_i_2d(ji,jl)  
     523               ht_i_2d(ji,jl)  =  v_i_2d(ji,jl) / a_i_2d(ji,jl)  
    524524               t_su_2d(ji,jl)  =  zaTsfn(ji,jl) / a_i_2d(ji,jl)  
    525525            ELSE 
    … …  
    529529         END DO 
    530530      END DO 
    531        
    532       CALL tab_2d_3d( nidx, idxice(1:nidx), ht_i_2d (1:nidx,1:jpl), ht_i   ) 
    533       CALL tab_2d_3d( nidx, idxice(1:nidx), a_i_2d  (1:nidx,1:jpl), a_i    ) 
    534       CALL tab_2d_3d( nidx, idxice(1:nidx), v_i_2d  (1:nidx,1:jpl), v_i    ) 
    535       CALL tab_2d_3d( nidx, idxice(1:nidx), v_s_2d  (1:nidx,1:jpl), v_s    ) 
    536       CALL tab_2d_3d( nidx, idxice(1:nidx), oa_i_2d (1:nidx,1:jpl), oa_i   ) 
    537       CALL tab_2d_3d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i  ) 
    538       CALL tab_2d_3d( nidx, idxice(1:nidx), a_ip_2d (1:nidx,1:jpl), a_ip   ) 
    539       CALL tab_2d_3d( nidx, idxice(1:nidx), v_ip_2d (1:nidx,1:jpl), v_ip   ) 
    540       CALL tab_2d_3d( nidx, idxice(1:nidx), t_su_2d (1:nidx,1:jpl), t_su   ) 
    541  
     531      ! 
     532      CALL tab_2d_3d( nidx, idxice(1:nidx), ht_i_2d (1:nidx,1:jpl), ht_i  ) 
     533      CALL tab_2d_3d( nidx, idxice(1:nidx), a_i_2d  (1:nidx,1:jpl), a_i   ) 
     534      CALL tab_2d_3d( nidx, idxice(1:nidx), v_i_2d  (1:nidx,1:jpl), v_i   ) 
     535      CALL tab_2d_3d( nidx, idxice(1:nidx), v_s_2d  (1:nidx,1:jpl), v_s   ) 
     536      CALL tab_2d_3d( nidx, idxice(1:nidx), oa_i_2d (1:nidx,1:jpl), oa_i  ) 
     537      CALL tab_2d_3d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i ) 
     538      CALL tab_2d_3d( nidx, idxice(1:nidx), a_ip_2d (1:nidx,1:jpl), a_ip  ) 
     539      CALL tab_2d_3d( nidx, idxice(1:nidx), v_ip_2d (1:nidx,1:jpl), v_ip  ) 
     540      CALL tab_2d_3d( nidx, idxice(1:nidx), t_su_2d (1:nidx,1:jpl), t_su  ) 
    542541      ! 
    543542   END SUBROUTINE ice_itd_shiftice 
    … …  
    559558      REAL(wp), DIMENSION(jpij,jpl-1) ::   zdaice, zdvice   ! ice area and volume transferred 
    560559      !!------------------------------------------------------------------ 
    561        
    562       DO jl = 1, jpl-1 
    563          jdonor(:,jl) = 0 
    564          zdaice(:,jl) = 0._wp 
    565          zdvice(:,jl) = 0._wp 
    566       END DO 
    567  
    568       !--------------------------------------- 
    569       ! identify thicknesses that are too big 
    570       !--------------------------------------- 
    571       DO jl = 1, jpl - 1 
    572  
    573          nidx = 0 ; idxice(:) = 0 
     560      ! 
     561      jdonor(:,:) = 0 
     562      zdaice(:,:) = 0._wp 
     563      zdvice(:,:) = 0._wp 
     564      ! 
     565      !                       !--------------------------------------- 
     566      DO jl = 1, jpl-1        ! identify thicknesses that are too big 
     567         !                    !--------------------------------------- 
     568         nidx = 0   ;   idxice(:) = 0 
    574569         DO jj = 1, jpj 
    575570            DO ji = 1, jpi 
    … …  
    578573                  idxice( nidx ) = (jj - 1) * jpi + ji                   
    579574               ENDIF 
    580             ENDDO 
    581          ENDDO 
    582           
     575            END DO 
     576         END DO 
     577         ! 
    583578!!clem   CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,jl)   ) 
    584579         CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl)    ) 
    585580         CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl)    ) 
    586  
     581         ! 
    587582         DO ji = 1, nidx 
    588583            jdonor(ji,jl)  = jl  
    … …  
    597592            zdaice(ji,jl)  = a_i_1d(ji) * 0.5_wp 
    598593            zdvice(ji,jl)  = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp 
    599              
    600          END DO 
    601  
     594         END DO 
     595         ! 
    602596         IF( nidx > 0 ) THEN 
    603597            CALL ice_itd_shiftice( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) )  ! Shift jl=>jl+1 
    … …  
    610604      END DO 
    611605 
    612       !----------------------------------------- 
    613       ! Identify thicknesses that are too small 
    614       !----------------------------------------- 
    615       DO jl = jpl - 1, 1, -1 
    616  
     606      !                       !----------------------------------------- 
     607      DO jl = jpl-1, 1, -1    ! Identify thicknesses that are too small 
     608         !                    !----------------------------------------- 
    617609         nidx = 0 ; idxice(:) = 0 
    618610         DO jj = 1, jpj 
    … …  
    622614                  idxice( nidx ) = (jj - 1) * jpi + ji                   
    623615               ENDIF 
    624             ENDDO 
    625          ENDDO 
    626  
     616            END DO 
     617         END DO 
     618         ! 
    627619         CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl+1)    ) ! jl+1 is ok 
    628620         CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl+1)    ) ! jl+1 is ok 
    … …  
    632624            zdvice(ji,jl) = v_i_1d(ji) 
    633625         END DO 
    634           
     626         ! 
    635627         IF( nidx > 0 ) THEN 
    636628            CALL ice_itd_shiftice( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) )  ! Shift jl+1=>jl 
    … …  
    640632            zdvice(1:nidx,jl) = 0._wp 
    641633         ENDIF 
    642  
     634         ! 
    643635      END DO 
    644636      ! 
    645637   END SUBROUTINE ice_itd_reb 
    646638 
     639#else 
     640   !!---------------------------------------------------------------------- 
     641   !!   Default option :         Empty module          NO LIM sea-ice model 
     642   !!---------------------------------------------------------------------- 
    647643#endif 
     644 
    648645   !!====================================================================== 
    649646END MODULE iceitd 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90

    r8426 r8486  
    1212   !!   'key_lim3'                                      LIM-3 sea-ice model 
    1313   !!---------------------------------------------------------------------- 
    14    USE par_oce          ! ocean parameters 
    15    USE dom_oce          ! ocean domain 
    16    USE phycst           ! physical constants (ocean directory)  
    17    USE sbc_oce, ONLY: sss_m, sst_m          ! surface boundary condition: ocean fields 
    18    USE ice1D            ! LIM thermodynamics 
    19    USE ice              ! LIM variables 
    20    USE icevar           ! LIM 
    21    USE icectl           ! control prints 
     14   USE par_oce        ! ocean parameters 
     15   USE dom_oce        ! ocean domain 
     16   USE phycst         ! physical constants (ocean directory)  
     17   USE sbc_oce , ONLY : sss_m, sst_m   ! surface boundary condition: ocean fields 
     18   USE ice1D          ! sea-ice: thermodynamics 
     19   USE ice            ! sea-ice: variables 
     20   USE icevar         ! sea-ice: operations 
     21   USE icectl         ! sea-ice: control prints 
    2222   ! 
    23    USE lbclnk           ! lateral boundary condition - MPP exchanges 
    24    USE lib_mpp          ! MPP library 
    25    USE in_out_manager   ! I/O manager 
    26    USE iom              ! I/O manager 
    27    USE lib_fortran      ! glob_sum 
    28    USE timing           ! Timing 
     23   USE lbclnk         ! lateral boundary condition - MPP exchanges 
     24   USE lib_mpp        ! MPP library 
     25   USE in_out_manager ! I/O manager 
     26   USE iom            ! I/O manager 
     27   USE lib_fortran    ! glob_sum 
     28   USE timing         ! Timing 
    2929 
    3030   IMPLICIT NONE 
    … …  
    3232 
    3333   PUBLIC   ice_rdgrft               ! called by ice_stp 
    34    PUBLIC   ice_rdgrft_icestrength 
    35    PUBLIC   ice_rdgrft_init 
    36    PUBLIC   ice_rdgrft_alloc        ! called by ice_init  
     34   PUBLIC   ice_rdgrft_icestrength   ! called by icerhg_evp 
     35   PUBLIC   ice_rdgrft_init          ! called by ice_stp 
     36   PUBLIC   ice_rdgrft_alloc         ! called by ice_init  
    3737 
    3838   !----------------------------------------------------------------------- 
    … …  
    5252   REAL(wp), PARAMETER ::   kraft   = 0.5_wp    ! rafting multipliyer 
    5353 
    54    REAL(wp) ::   Cp                             !  
     54!!gm Cp is  1) not DOCTOR,  
     55!!          2) misleading name: heat capacity instead of a constant, 
     56!!          3) recomputed at each time-step, whereas it is stored in the module memory... 
     57!!             ===>>> compute it one for all inside the IF( kt == nit000 ) (i.e. without the ".AND. lwp") 
     58   REAL(wp)            ::   Cp                  ! ??? !!gm  Not doctor !  
     59    
    5560   ! 
    5661   ! 
    5762   !!---------------------------------------------------------------------- 
    58    !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010) 
     63   !! NEMO/ICE 4.0 , NEMO Consortium (2017) 
    5964   !! $Id: icerdgrft.F90 8378 2017-07-26 13:55:59Z clem $ 
    6065   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    6671      !!                ***  ROUTINE ice_rdgrft_alloc *** 
    6772      !!---------------------------------------------------------------------! 
    68       ALLOCATE(                                                                      & 
    69          !* Variables shared among ridging subroutines 
    70          &      asum (jpi,jpj)     , athorn(jpi,jpj,0:jpl) , aksum (jpi,jpj)     ,   & 
    71          &      hrmin(jpi,jpj,jpl) , hraft(jpi,jpj,jpl)    , aridge(jpi,jpj,jpl) ,   & 
    72          &      hrmax(jpi,jpj,jpl) , krdg (jpi,jpj,jpl)    , araft (jpi,jpj,jpl) , STAT=ice_rdgrft_alloc ) 
     73      ALLOCATE( asum (jpi,jpj)     , athorn(jpi,jpj,0:jpl) , aksum (jpi,jpj)     ,     & 
     74         &      hrmin(jpi,jpj,jpl) , hraft (jpi,jpj,jpl)   , aridge(jpi,jpj,jpl) ,     & 
     75         &      hrmax(jpi,jpj,jpl) , krdg  (jpi,jpj,jpl)   , araft (jpi,jpj,jpl) , STAT=ice_rdgrft_alloc ) 
    7376         ! 
    7477      IF( ice_rdgrft_alloc /= 0 )   CALL ctl_warn( 'ice_rdgrft_alloc: failed to allocate arrays' ) 
    … …  
    105108      INTEGER, INTENT(in) ::   kt     ! number of iteration 
    106109      !! 
    107       INTEGER  ::   ji, jj, jk, jl        ! dummy loop index 
    108       INTEGER  ::   niter                 ! local integer  
    109       INTEGER  ::   iterate_ridging       ! if true, repeat the ridging 
    110       REAL(wp) ::   za, zfac              ! local scalar 
     110      INTEGER  ::   ji, jj, jk, jl     ! dummy loop index 
     111      INTEGER  ::   niter              ! local integer  
     112      INTEGER  ::   iterate_ridging    ! if =1, repeat the ridging 
     113      REAL(wp) ::   za, zfac, zcs_2    ! local scalar 
    111114      CHARACTER (len = 15) ::   fieldid 
    112       REAL(wp), DIMENSION(jpi,jpj)   ::   closing_net     ! net rate at which area is removed    (1/s) 
    113                                                                ! (ridging ice area - area of new ridges) / dt 
    114       REAL(wp), DIMENSION(jpi,jpj)   ::   divu_adv        ! divu as implied by transport scheme  (1/s) 
    115       REAL(wp), DIMENSION(jpi,jpj)   ::   opning          ! rate of opening due to divergence/shear 
    116       REAL(wp), DIMENSION(jpi,jpj)   ::   closing_gross   ! rate at which area removed, not counting area of new ridges 
     115      REAL(wp), DIMENSION(jpi,jpj) ::   closing_net     ! net rate at which area is removed    (1/s) 
     116      !                                                 ! (ridging ice area - area of new ridges) / dt 
     117      REAL(wp), DIMENSION(jpi,jpj) ::   divu_adv        ! divu as implied by transport scheme  (1/s) 
     118      REAL(wp), DIMENSION(jpi,jpj) ::   opning          ! rate of opening due to divergence/shear 
     119      REAL(wp), DIMENSION(jpi,jpj) ::   closing_gross   ! rate at which area removed, not counting area of new ridges 
    117120      ! 
    118121      INTEGER, PARAMETER ::   nitermax = 20     
    … …  
    124127      IF( kt == nit000 .AND. lwp ) THEN 
    125128         WRITE(numout,*) 
    126          WRITE(numout,*)'icerdgrft' 
    127          WRITE(numout,*)'~~~~~~~~~' 
     129         WRITE(numout,*)'icerdgrft : ice ridging and rafting' 
     130         WRITE(numout,*)'~~~~~~~~~~' 
    128131      ENDIF 
    129  
    130       ! conservation test 
    131       IF( ln_limdiachk ) CALL ice_cons_hsm(0, 'icerdgrft', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     132!!gm should be:       
     133!      IF( kt == nit000 ) THEN 
     134!         IF(lwp) WRITE(numout,*) 
     135!         IF(lwp) WRITE(numout,*)'icerdgrft : ???' 
     136!         IF(lwp) WRITE(numout,*)'~~~~~~~~~~' 
     137!         ! 
     138!         Cp    = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0      ! proport const for PE 
     139!         ! 
     140!!gm why not adding here zcs_2 computation 
     141!         ! 
     142!      ENDIF 
     143!!gm end 
     144       
     145      !                    ! conservation test 
     146      IF( ln_limdiachk )   CALL ice_cons_hsm(0, 'icerdgrft', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    132147 
    133148      !-----------------------------------------------------------------------------! 
    134149      ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons 
    135150      !-----------------------------------------------------------------------------! 
    136       Cp = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0             ! proport const for PE 
    137       ! 
    138       CALL ice_rdgrft_ridgeprep                                    ! prepare ridging 
    139       ! 
    140  
    141       DO jj = 1, jpj                                     ! Initialize arrays. 
     151      Cp    = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0      ! proport const for PE 
     152      zcs_2 = rn_cs * 0.5_wp 
     153      ! 
     154      CALL ice_rdgrft_ridgeprep                             ! prepare ridging 
     155      ! 
     156      DO jj = 1, jpj                                        ! Initialize arrays. 
    142157         DO ji = 1, jpi 
    143  
    144158            !-----------------------------------------------------------------------------! 
    145159            ! 2) Dynamical inputs (closing rate, divu_adv, opning) 
    … …  
    161175            !  (thick, newly ridged ice). 
    162176 
    163             closing_net(ji,jj) = rn_cs * 0.5 * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp ) 
     177            closing_net(ji,jj) = zcs_2 * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp ) 
    164178 
    165179            ! 2.2 divu_adv 
    … …  
    233247         ! 3.3 Redistribute area, volume, and energy. 
    234248         !-----------------------------------------------------------------------------! 
    235  
    236249         CALL ice_rdgrft_ridgeshift( opning, closing_gross ) 
    237  
    238250          
    239251         ! 3.4 Compute total area of ice plus open water after ridging. 
    … …  
    246258         ! Check whether asum = 1.  If not (because the closing and opening 
    247259         ! rates were reduced above), ridge again with new rates. 
    248  
    249260         iterate_ridging = 0 
    250261         DO jj = 1, jpj 
    … …  
    262273            END DO 
    263274         END DO 
    264  
    265275         IF( lk_mpp )   CALL mpp_max( iterate_ridging ) 
    266276 
    267277         ! Repeat if necessary. 
    268278         ! NOTE: If strength smoothing is turned on, the ridging must be 
    269          ! iterated globally because of the boundary update in the  
    270          ! smoothing. 
    271  
     279         ! iterated globally because of the boundary update in the smoothing. 
    272280         niter = niter + 1 
    273  
     281         ! 
    274282         IF( iterate_ridging == 1 ) THEN 
    275283            CALL ice_rdgrft_ridgeprep 
    … …  
    279287            ENDIF 
    280288         ENDIF 
    281  
     289         ! 
    282290      END DO !! on the do while over iter 
    283291 
    … …  
    287295      ! control prints 
    288296      !-----------------------------------------------------------------------------! 
    289       ! conservation test 
    290       IF( ln_limdiachk ) 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' ) 
     297      !                    ! conservation test 
     298      IF( ln_limdiachk )   CALL ice_cons_hsm(1, 'icerdgrft', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     299 
     300      !                    ! control prints 
     301      IF( ln_ctl       )   CALL ice_prt3D( 'icerdgrft' ) 
    294302      ! 
    295303      IF( nn_timing == 1 )  CALL timing_stop('icerdgrft') 
     304      ! 
    296305   END SUBROUTINE ice_rdgrft 
     306 
    297307 
    298308   SUBROUTINE ice_rdgrft_ridgeprep 
    … …  
    305315      !!              participating in ridging and of the resulting ridges. 
    306316      !!---------------------------------------------------------------------! 
    307       INTEGER  ::   ji,jj, jl    ! dummy loop indices 
    308       REAL(wp) ::   Gstari, astari, hrmean, zdummy   ! local scalar 
     317      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
     318      REAL(wp) ::   Gstari, astari, hrmean, zdummy   ! local scalar     !!gm DOCTOR norme should start with z !!!!! 
    309319      REAL(wp), DIMENSION(jpi,jpj,-1:jpl) ::   Gsum      ! Gsum(n) = sum of areas in categories 0 to n 
    310320      !------------------------------------------------------------------------------! 
    311321 
    312       Gstari     = 1.0/rn_gstar     
    313       astari     = 1.0/rn_astar     
    314       aksum(:,:)    = 0.0 
    315       athorn(:,:,:) = 0.0 
    316       aridge(:,:,:) = 0.0 
    317       araft (:,:,:) = 0.0 
    318  
    319       ! Zero out categories with very small areas 
    320       CALL ice_var_zapsmall 
    321  
    322       ! Ice thickness needed for rafting 
    323       DO jl = 1, jpl 
    324          DO jj = 1, jpj 
    325             DO ji = 1, jpi 
    326                rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) 
    327                ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 
    328            END DO 
    329          END DO 
    330       END DO 
     322      Gstari        = 1._wp / rn_gstar     
     323      astari        = 1._wp / rn_astar     
     324      aksum(:,:)    = 0._wp 
     325      athorn(:,:,:) = 0._wp 
     326      aridge(:,:,:) = 0._wp 
     327      araft (:,:,:) = 0._wp 
     328 
     329      CALL ice_var_zapsmall   ! Zero out categories with very small areas 
     330 
     331!      DO jl = 1, jpl          ! Ice thickness needed for rafting 
     332!         DO jj = 1, jpj 
     333!            DO ji = 1, jpi 
     334!!gm               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) 
     335!!gm               ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 
     336!               IF( a_i(ji,jj,jl) >= epsi20 ) THEN   ;   ht_i(ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl) 
     337!               ELSE                                 ;   ht_i(ji,jj,jl) = 0._wp 
     338!               ENDIF 
     339!           END DO 
     340!         END DO 
     341!      END DO 
     342!!gm or even better : 
     343!     !                       ! Ice thickness needed for rafting 
     344      WHERE( a_i(:,:,:) >= epsi20 )   ;   ht_i(:,:,:) = v_i (:,:,:) / a_i(:,:,:) 
     345      ELSEWHERE                       ;   ht_i(:,:,:) = 0._wp 
     346      END WHERE 
     347!!gm end 
    331348 
    332349      !------------------------------------------------------------------------------! 
    333350      ! 1) Participation function  
    334351      !------------------------------------------------------------------------------! 
    335  
     352      ! 
    336353      ! Compute total area of ice plus open water. 
    337354      ! This is in general not equal to one because of divergence during transport 
    338355      asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 ) 
    339  
     356      ! 
    340357      ! Compute cumulative thickness distribution function 
    341358      ! Compute the cumulative thickness distribution function Gsum, 
    … …  
    348365         Gsum(:,:,jl) = Gsum(:,:,jl-1) + a_i(:,:,jl) 
    349366      END DO 
    350  
     367      ! 
    351368      ! Normalize the cumulative distribution to 1 
    352369      DO jl = 0, jpl 
    … …  
    366383      ! athorn is always >= 0 and SUM(athorn(0:jpl))=1 
    367384      !----------------------------------------------------------------- 
    368  
     385      ! 
    369386      IF( nn_partfun == 0 ) THEN       !--- Linear formulation (Thorndike et al., 1975) 
    370387         DO jl = 0, jpl     
    … …  
    383400            END DO 
    384401         END DO 
    385  
     402         ! 
    386403      ELSE                             !--- Exponential, more stable formulation (Lipscomb et al, 2007) 
    387404         !                         
    … …  
    396413      ENDIF 
    397414 
    398       ! --- Ridging and rafting participation concentrations --- ! 
    399       IF( ln_rafting .AND. ln_ridging ) THEN 
    400          ! 
     415      !                                !--- Ridging and rafting participation concentrations 
     416      IF( ln_rafting .AND. ln_ridging ) THEN             !- ridging & rafting 
    401417         DO jl = 1, jpl 
    402418            DO jj = 1, jpj  
    … …  
    408424            END DO 
    409425         END DO 
    410          ! 
    411       ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN 
    412          ! 
     426      ELSEIF( ln_ridging .AND. .NOT.ln_rafting ) THEN   !- ridging alone 
    413427         DO jl = 1, jpl 
    414428            aridge(:,:,jl) = athorn(:,:,jl) 
    415429         END DO 
    416          ! 
    417       ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN 
    418          ! 
     430      ELSEIF( ln_rafting .AND. .NOT.ln_ridging ) THEN   !- rafting alone    
    419431         DO jl = 1, jpl 
    420432            araft(:,:,jl) = athorn(:,:,jl) 
    421433         END DO 
    422          ! 
    423434      ENDIF 
    424435 
    … …  
    454465         DO jj = 1, jpj 
    455466            DO ji = 1, jpi 
    456                 
    457                IF( athorn(ji,jj,jl) > 0._wp ) THEN 
     467               IF ( athorn(ji,jj,jl) > 0._wp ) THEN 
    458468                  hrmean          = MAX( SQRT( rn_hstar * ht_i(ji,jj,jl) ), ht_i(ji,jj,jl) * krdgmin ) 
    459469                  hrmin(ji,jj,jl) = MIN( 2._wp * ht_i(ji,jj,jl), 0.5_wp * ( hrmean + ht_i(ji,jj,jl) ) ) 
    460470                  hrmax(ji,jj,jl) = 2._wp * hrmean - hrmin(ji,jj,jl) 
    461471                  hraft(ji,jj,jl) = ht_i(ji,jj,jl) / kraft 
    462                   krdg(ji,jj,jl) = ht_i(ji,jj,jl) / MAX( hrmean, epsi20 ) 
    463  
     472                  krdg (ji,jj,jl) = ht_i(ji,jj,jl) / MAX( hrmean, epsi20 ) 
     473                  ! 
    464474                  ! Normalization factor : aksum, ensures mass conservation 
    465475                  aksum(ji,jj) = aksum(ji,jj) + aridge(ji,jj,jl) * ( 1._wp - krdg(ji,jj,jl) )    & 
    466476                     &                        + araft (ji,jj,jl) * ( 1._wp - kraft          ) 
    467  
    468477               ELSE 
    469478                  hrmin(ji,jj,jl)  = 0._wp  
    … …  
    472481                  krdg (ji,jj,jl)  = 1._wp 
    473482               ENDIF 
    474  
    475483            END DO 
    476484         END DO 
    477485      END DO 
    478       ! 
    479486      ! 
    480487   END SUBROUTINE ice_rdgrft_ridgeprep 
    … …  
    493500      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   closing_gross  ! rate at which area removed, excluding area of new ridges 
    494501      ! 
    495       CHARACTER (len=80) ::   fieldid   ! field identifier 
     502      CHARACTER (len=80) ::   fieldid   ! field identifier     !!gm DOCTOR name wrong 
    496503      ! 
    497504      INTEGER ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices 
    … …  
    740747 
    741748               ! Compute the fraction of rafted ice area and volume going to thickness category jl2 
    742                IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) >  hi_max(jl2-1) ) THEN 
    743                   zswitch(ij) = 1._wp 
    744                ELSE 
    745                   zswitch(ij) = 0._wp                   
     749!!gm see above               IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) >  hi_max(jl2-1) ) THEN 
     750               IF( hi_max(jl2-1) < hraft(ji,jj,jl1) .AND. hraft(ji,jj,jl1) <= hi_max(jl2)  ) THEN   ;   zswitch(ij) = 1._wp 
     751               ELSE                                                                                 ;   zswitch(ij) = 0._wp                   
    746752               ENDIF 
    747  
     753               ! 
    748754               a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + ( ardg2 (ij) * farea    + arft2 (ij) * zswitch(ij) ) 
    749755               oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + ( oirdg2(ij) * farea    + oirft2(ij) * zswitch(ij) ) 
    … …  
    756762               ! MV MP 2016 
    757763               IF ( nn_pnd_scheme > 0 ) THEN 
    758                   v_ip (ji,jj,jl2) = v_ip (ji,jj,jl2)  + ( vprdg (ij) * rn_fpondrdg * fvol(ij)  +   & 
    759                                                        &   vprft (ij) * rn_fpondrft * zswitch(ij) ) 
    760                   a_ip (ji,jj,jl2) = a_ip(ji,jj,jl2)   + ( aprdg2(ij) * rn_fpondrdg * farea +       &  
    761                                                        &   aprft2(ij) * rn_fpondrft * zswitch(ji) ) 
     764                  v_ip (ji,jj,jl2) = v_ip(ji,jj,jl2) + (   vprdg (ij) * rn_fpondrdg * fvol   (ij)   & 
     765                     &                                   + vprft (ij) * rn_fpondrft * zswitch(ij)  ) 
     766                  a_ip (ji,jj,jl2) = a_ip(ji,jj,jl2) + (   aprdg2(ij) * rn_fpondrdg * farea         &  
     767                     &                                   + aprft2(ij) * rn_fpondrft * zswitch(ji)  ) 
    762768               ENDIF 
    763769               ! END MV MP 2016 
    764  
    765770            END DO 
    766771 
    … …  
    774779            ! 
    775780         END DO ! jl2 
    776           
     781         ! 
    777782      END DO ! jl1 (deforming categories) 
    778783 
    … …  
    782787      ! 
    783788   END SUBROUTINE ice_rdgrft_ridgeshift 
     789 
    784790 
    785791   SUBROUTINE ice_rdgrft_icestrength( kstrngth ) 
    … …  
    798804      !!---------------------------------------------------------------------- 
    799805      INTEGER, INTENT(in) ::   kstrngth    ! = 1 for Rothrock formulation, 0 for Hibler (1979) 
     806      ! 
    800807      INTEGER             ::   ji,jj, jl   ! dummy loop indices 
    801       INTEGER             ::   ksmooth     ! smoothing the resistance to deformation 
    802       INTEGER             ::   numts_rm    ! number of time steps for the P smoothing 
     808      INTEGER             ::   ksmooth     ! smoothing the resistance to deformation    !!gm not DOCTOR : start with i !!! 
     809      INTEGER             ::   numts_rm    ! number of time steps for the P smoothing    !!gm not DOCTOR : start with i !!! 
    803810      REAL(wp)            ::   zp, z1_3    ! local scalars 
    804811      REAL(wp), DIMENSION(jpi,jpj) ::   zworka           ! temporary array used here 
    … …  
    880887      ! 6) Smoothing ice strength 
    881888      !------------------------------------------------------------------------------! 
    882       ! 
    883       !------------------- 
    884       ! Spatial smoothing 
    885       !------------------- 
    886       IF ( ksmooth == 1 ) THEN 
    887  
     889      SELECT CASE( ksmooth ) 
     890      !                       !------------------- 
     891      CASE( 1 )               ! Spatial smoothing 
     892         !                    !------------------- 
    888893         DO jj = 2, jpjm1 
    889894            DO ji = 2, jpim1 
    … …  
    905910         END DO 
    906911         CALL lbc_lnk( strength, 'T', 1. ) 
    907  
    908       ENDIF 
    909  
    910       !-------------------- 
    911       ! Temporal smoothing 
    912       !-------------------- 
    913       IF ( ksmooth == 2 ) THEN 
    914  
     912         ! 
     913         !                    !-------------------- 
     914      CASE( 2 )               ! Temporal smoothing 
     915         !                    !-------------------- 
    915916         IF ( kt_ice == nit000 ) THEN 
    916917            zstrp1(:,:) = 0._wp 
    917918            zstrp2(:,:) = 0._wp 
    918919         ENDIF 
    919  
     920         ! 
    920921         DO jj = 2, jpjm1 
    921922            DO ji = 2, jpim1 
    … …  
    925926                  IF ( zstrp2(ji,jj) > 0._wp ) numts_rm = numts_rm + 1 
    926927                  zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / numts_rm 
    927                   zstrp2(ji,jj) = zstrp1(ji,jj) 
    928                   zstrp1(ji,jj) = strength(ji,jj) 
     928                  zstrp2  (ji,jj) = zstrp1  (ji,jj) 
     929                  zstrp1  (ji,jj) = strength(ji,jj) 
    929930                  strength(ji,jj) = zp 
    930931               ENDIF 
    931932            END DO 
    932933         END DO 
    933  
    934934         CALL lbc_lnk( strength, 'T', 1. )      ! Boundary conditions 
    935  
    936       ENDIF ! ksmooth 
     935         ! 
     936      END SELECT 
    937937      ! 
    938938   END SUBROUTINE ice_rdgrft_icestrength 
     939 
    939940 
    940941   SUBROUTINE ice_rdgrft_init 
    941942      !!------------------------------------------------------------------- 
    942       !!                   ***  ROUTINE ice_rdgrft_init *** 
     943      !!                  ***  ROUTINE ice_rdgrft_init *** 
    943944      !! 
    944945      !! ** Purpose :   Physical constants and parameters linked  
    … …  
    952953      !!------------------------------------------------------------------- 
    953954      INTEGER :: ios                 ! Local integer output status for namelist read 
    954       NAMELIST/namiceitdme/ rn_cs, nn_partfun, rn_gstar, rn_astar,             &  
    955         &                   ln_ridging, rn_hstar, rn_por_rdg, rn_fsnowrdg, rn_fpondrdg, &  
    956                             ln_rafting, rn_hraft, rn_craft,   rn_fsnowrft, rn_fpondrft 
     955      !! 
     956      NAMELIST/namiceitdme/ rn_cs     , nn_partfun, rn_gstar  , rn_astar   ,                &  
     957         &                  ln_ridging, rn_hstar  , rn_por_rdg, rn_fsnowrdg, rn_fpondrdg,   &  
     958         &                  ln_rafting, rn_hraft  , rn_craft  , rn_fsnowrft, rn_fpondrft 
    957959      !!------------------------------------------------------------------- 
    958960      ! 
    … …  
    960962      READ  ( numnam_ice_ref, namiceitdme, IOSTAT = ios, ERR = 901) 
    961963901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in reference namelist', lwp ) 
    962  
     964      ! 
    963965      REWIND( numnam_ice_cfg )              ! Namelist namiceitdme in configuration namelist : Ice mechanical ice redistribution 
    964966      READ  ( numnam_ice_cfg, namiceitdme, IOSTAT = ios, ERR = 902 ) 
    … …  
    992994   !!   Default option         Empty module          NO LIM-3 sea-ice model 
    993995   !!---------------------------------------------------------------------- 
    994 CONTAINS 
    995    SUBROUTINE ice_rdgrft           ! Empty routines 
    996    END SUBROUTINE ice_rdgrft 
    997    SUBROUTINE ice_rdgrft_icestrength 
    998    END SUBROUTINE ice_rdgrft_icestrength 
    999    SUBROUTINE ice_rdgrft_init 
    1000    END SUBROUTINE ice_rdgrft_init 
    1001996#endif 
     997 
    1002998   !!====================================================================== 
    1003999END MODULE icerdgrft 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg.F90

    r8426 r8486  
    1010#if defined key_lim3 
    1111   !!---------------------------------------------------------------------- 
    12    !!   'key_lim3' :                                 LIM3 sea-ice model 
     12   !!   'key_lim3'                                       LIM3 sea-ice model 
    1313   !!---------------------------------------------------------------------- 
    1414   !!    ice_rhg      : computes ice velocities 
    1515   !!    ice_rhg_init : initialization and namelist read 
    1616   !!---------------------------------------------------------------------- 
    17    USE phycst           ! physical constants 
    18    USE dom_oce          ! ocean space and time domain 
    19    USE ice              ! LIM-3 variables 
    20    USE icerhg_evp       ! EVP rheology 
    21    USE icectl           ! control prints 
    22    USE icevar 
     17   USE phycst         ! physical constants 
     18   USE dom_oce        ! ocean space and time domain 
     19   USE ice            ! sea-ice: variables 
     20   USE icerhg_evp     ! sea-ice: EVP rheology 
     21   USE icectl         ! sea-ice: control prints 
     22   USE icevar         ! sea-ice: operations 
    2323   ! 
    24    USE lbclnk           ! lateral boundary conditions - MPP exchanges 
    25    USE lib_mpp          ! MPP library 
    26    USE in_out_manager   ! I/O manager 
    27    USE lib_fortran      ! glob_sum 
    28    USE timing           ! Timing 
     24   USE lbclnk         ! lateral boundary conditions - MPP exchanges 
     25   USE lib_mpp        ! MPP library 
     26   USE in_out_manager ! I/O manager 
     27   USE lib_fortran    ! glob_sum 
     28   USE timing         ! Timing 
    2929 
    3030   IMPLICIT NONE 
    … …  
    3737#  include "vectopt_loop_substitute.h90" 
    3838   !!---------------------------------------------------------------------- 
    39    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     39   !! NEMO/ICE 4.0 , NEMO Consortium (2017) 
    4040   !! $Id: icerhg.F90 8378 2017-07-26 13:55:59Z clem $ 
    4141   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    5353      !!                      - shear, divergence and delta (shear_i, divu_i, delta_i) 
    5454      !!-------------------------------------------------------------------- 
    55       INTEGER, INTENT(in) ::   kt     ! number of iteration 
     55      INTEGER, INTENT(in) ::   kt     ! ice time step 
    5656      !! 
    57       INTEGER  :: jl ! dummy loop indices 
     57      INTEGER  ::   jl  ! dummy loop indices 
    5858      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    5959      !!-------------------------------------------------------------------- 
    60  
     60      ! 
    6161      IF( nn_timing == 1 )  CALL timing_start('icerhg') 
    62  
     62      ! 
    6363      IF( kt == nit000 .AND. lwp ) THEN 
    6464         WRITE(numout,*) 
    65          WRITE(numout,*)'icerhg' 
    66          WRITE(numout,*)'~~~~~~' 
     65         WRITE(numout,*)'ice_rhg : sea-ice rheology' 
     66         WRITE(numout,*)'~~~~~~~~' 
    6767      ENDIF 
    6868 
    69       CALL ice_var_agg(1)   ! aggregate ice categories 
     69      CALL ice_var_agg(1)           ! -- aggregate ice categories 
    7070      ! 
    71       ! conservation test 
    72       IF( ln_limdiachk ) CALL ice_cons_hsm(0, 'icerhg', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    73        
    74       ! Landfast ice parameterization: define max bottom friction 
    75       tau_icebfr(:,:) = 0._wp 
    76       IF( ln_landfast ) THEN 
     71      !                             ! -- conservation test 
     72      IF( ln_limdiachk )   CALL ice_cons_hsm(0, 'icerhg', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     73      !                       
     74      IF( ln_landfast ) THEN        ! -- Landfast ice parameterization: define max bottom friction 
    7775         DO jl = 1, jpl 
    78             WHERE( ht_i(:,:,jl) > ht_n(:,:) * rn_gamma )  tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 
     76            WHERE( ht_i(:,:,jl) > ht_n(:,:) * rn_gamma )   ;   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 
     77            ELSEWHERE                                      ;   tau_icebfr(:,:) = 0._wp 
     78            END WHERE 
    7979         END DO 
    8080      ENDIF 
    … …  
    8383      ! Rheology (ice dynamics) 
    8484      ! -----------------------    
    85       IF( nn_limdyn /= 0 ) THEN                          ! -- Ice dynamics 
    86  
     85      IF( nn_limdyn /= 0 ) THEN     ! -- Ice dynamics 
     86         ! 
    8787         CALL ice_rhg_evp( stress1_i, stress2_i, stress12_i, u_ice, v_ice, shear_i, divu_i, delta_i ) 
    88  
    89       ELSE 
    90  
    91          u_ice(:,:) = rn_uice * umask(:,:,1)             !     or prescribed velocity 
     88         ! 
     89      ELSE                          ! -- prescribed uniform velocity 
     90         ! 
     91         u_ice(:,:) = rn_uice * umask(:,:,1) 
    9292         v_ice(:,:) = rn_vice * vmask(:,:,1) 
    9393         !!CALL RANDOM_NUMBER(u_ice(:,:)) 
    9494         !!CALL RANDOM_NUMBER(v_ice(:,:)) 
    95  
     95         ! 
    9696      ENDIF 
    9797      ! 
    98       ! conservation test 
    99       IF( ln_limdiachk ) CALL ice_cons_hsm(1, 'icerhg', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    100  
    101       ! Control prints 
    102       IF( ln_ctl )       CALL ice_prt3D( 'icerhg' ) 
     98      !                                                   !- conservation test 
     99      IF( ln_limdiachk   )   CALL ice_cons_hsm(1, 'icerhg', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     100      IF( ln_ctl         )   CALL ice_prt3D  ('icerhg')   !- Control prints 
     101      IF( nn_timing == 1 )   CALL timing_stop('icerhg')   !- timing 
    103102      ! 
    104       IF( nn_timing == 1 )  CALL timing_stop('icerhg') 
    105  
    106103   END SUBROUTINE ice_rhg 
    107104 
    … …  
    119116      !! ** input   :   Namelist namicedyn 
    120117      !!------------------------------------------------------------------- 
    121       INTEGER  ::   ios                 ! Local integer output status for namelist read 
    122       NAMELIST/namicedyn/ nn_limadv, nn_limadv_ord,                                & 
    123          &                nn_icestr, rn_pe_rdg, rn_pstar, rn_crhg, ln_icestr_bvf,  & 
    124          &                rn_ishlat, rn_cio, rn_creepl, rn_ecc,                    & 
    125          &                nn_nevp, rn_relast, ln_landfast, rn_gamma, rn_icebfr, rn_lfrelax 
     118      INTEGER ::   ios   ! Local integer output status for namelist read 
     119      !! 
     120      NAMELIST/namicedyn/ nn_limadv  , nn_limadv_ord,                                       & 
     121         &                nn_icestr  , rn_pe_rdg, rn_pstar , rn_crhg, ln_icestr_bvf     ,   & 
     122         &                rn_ishlat  , rn_cio   , rn_creepl, rn_ecc , nn_nevp, rn_relast,   & 
     123         &                ln_landfast, rn_gamma , rn_icebfr, rn_lfrelax 
    126124      !!------------------------------------------------------------------- 
    127  
    128       REWIND( numnam_ice_ref )              ! Namelist namicedyn in reference namelist : Ice dynamics 
     125      ! 
     126      REWIND( numnam_ice_ref )         ! Namelist namicedyn in reference namelist : Ice dynamics 
    129127      READ  ( numnam_ice_ref, namicedyn, IOSTAT = ios, ERR = 901) 
    130128901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicedyn in reference namelist', lwp ) 
    131  
    132       REWIND( numnam_ice_cfg )              ! Namelist namicedyn in configuration namelist : Ice dynamics 
     129      ! 
     130      REWIND( numnam_ice_cfg )         ! Namelist namicedyn in configuration namelist : Ice dynamics 
    133131      READ  ( numnam_ice_cfg, namicedyn, IOSTAT = ios, ERR = 902 ) 
    134132902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicedyn in configuration namelist', lwp ) 
    135133      IF(lwm) WRITE ( numoni, namicedyn ) 
    136        
    137       IF(lwp) THEN                        ! control print 
     134      ! 
     135      IF(lwp) THEN                     ! control print 
    138136         WRITE(numout,*) 
    139137         WRITE(numout,*) 'ice_rhg_init : ice parameters for ice dynamics ' 
    140138         WRITE(numout,*) '~~~~~~~~~~~~' 
    141          ! limtrp 
    142          WRITE(numout,*)'    choose the advection scheme (-1=Prather, 0=Ulimate-Macho)   nn_limadv     = ', nn_limadv  
    143          WRITE(numout,*)'    choose the order of the scheme (if ultimate)                nn_limadv_ord = ', nn_limadv_ord   
    144          ! icerdgrft 
    145          WRITE(numout,*)'    ice strength parameterization (0=Hibler 1=Rothrock)         nn_icestr     = ', nn_icestr  
    146          WRITE(numout,*)'    Ratio of ridging work to PotEner change in ridging          rn_pe_rdg     = ', rn_pe_rdg  
    147          WRITE(numout,*) '   first bulk-rheology parameter                               rn_pstar      = ', rn_pstar 
    148          WRITE(numout,*) '   second bulk-rhelogy parameter                               rn_crhg       = ', rn_crhg 
    149          WRITE(numout,*)'    Including brine volume in ice strength comp.                ln_icestr_bvf = ', ln_icestr_bvf 
    150          ! icerhg_evp 
    151          WRITE(numout,*) '   lateral boundary condition for sea ice dynamics             rn_ishlat     = ', rn_ishlat 
    152          WRITE(numout,*) '   drag coefficient for oceanic stress                         rn_cio        = ', rn_cio 
    153          WRITE(numout,*) '   creep limit                                                 rn_creepl     = ', rn_creepl 
    154          WRITE(numout,*) '   eccentricity of the elliptical yield curve                  rn_ecc        = ', rn_ecc 
    155          WRITE(numout,*) '   number of iterations for subcycling                         nn_nevp       = ', nn_nevp 
    156          WRITE(numout,*) '   ratio of elastic timescale over ice time step               rn_relast     = ', rn_relast 
    157          WRITE(numout,*) '   Landfast: param (T or F)                                    ln_landfast   = ', ln_landfast 
    158          WRITE(numout,*) '      T: fraction of ocean depth that ice must reach           rn_gamma      = ', rn_gamma 
    159          WRITE(numout,*) '      T: maximum bottom stress per unit area of contact        rn_icebfr     = ', rn_icebfr 
    160          WRITE(numout,*) '      T: relax time scale (s-1) to reach static friction       rn_lfrelax    = ', rn_lfrelax 
     139         WRITE(numout,*) '   Namelist namicedyn' 
     140         WRITE(numout,*) '      advection scheme for ice transport (limtrp)' 
     141         WRITE(numout,*) '         type of advection scheme (-1=Prather, 0=Ulimate-Macho)   nn_limadv     = ', nn_limadv  
     142         WRITE(numout,*) '         order of the scheme for Ultimate-Macho case              nn_limadv_ord = ', nn_limadv_ord 
     143         WRITE(numout,*) '      ridging/rafting (icerdgrft)' 
     144         WRITE(numout,*) '         ice strength parameterization (0=Hibler 1=Rothrock)      nn_icestr     = ', nn_icestr  
     145         WRITE(numout,*) '         Ratio of ridging work to PotEner change in ridging       rn_pe_rdg     = ', rn_pe_rdg  
     146         WRITE(numout,*) '         1st bulk-rheology parameter                              rn_pstar      = ', rn_pstar 
     147         WRITE(numout,*) '         2nd bulk-rhelogy parameter                               rn_crhg       = ', rn_crhg 
     148         WRITE(numout,*) '         brine volume included in ice strength computation        ln_icestr_bvf = ', ln_icestr_bvf 
     149         WRITE(numout,*) '      rheology EVP (icerhg_evp)' 
     150         WRITE(numout,*) '         lateral boundary condition for sea ice dynamics          rn_ishlat     = ', rn_ishlat 
     151         WRITE(numout,*) '         drag coefficient for oceanic stress                      rn_cio        = ', rn_cio 
     152         WRITE(numout,*) '         creep limit                                              rn_creepl     = ', rn_creepl 
     153         WRITE(numout,*) '         eccentricity of the elliptical yield curve               rn_ecc        = ', rn_ecc 
     154         WRITE(numout,*) '         number of iterations for subcycling                      nn_nevp       = ', nn_nevp 
     155         WRITE(numout,*) '         ratio of elastic timescale over ice time step            rn_relast     = ', rn_relast 
     156         WRITE(numout,*) '      Landfast: param (T or F)                                    ln_landfast   = ', ln_landfast 
     157         WRITE(numout,*) '         fraction of ocean depth that ice must reach              rn_gamma      = ', rn_gamma 
     158         WRITE(numout,*) '         maximum bottom stress per unit area of contact           rn_icebfr     = ', rn_icebfr 
     159         WRITE(numout,*) '         relax time scale (s-1) to reach static friction          rn_lfrelax    = ', rn_lfrelax 
    161160      ENDIF 
    162161      ! 
    163       IF     (      rn_ishlat == 0.                ) THEN   ;   IF(lwp) WRITE(numout,*) '   ice lateral  free-slip ' 
    164       ELSEIF (      rn_ishlat == 2.                ) THEN   ;   IF(lwp) WRITE(numout,*) '   ice lateral  no-slip ' 
    165       ELSEIF ( 0. < rn_ishlat .AND. rn_ishlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ice lateral  partial-slip ' 
    166       ELSEIF ( 2. < rn_ishlat                      ) THEN   ;   IF(lwp) WRITE(numout,*) '   ice lateral  strong-slip ' 
     162      IF     (      rn_ishlat == 0.                ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  free-slip' 
     163      ELSEIF (      rn_ishlat == 2.                ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  no-slip' 
     164      ELSEIF ( 0. < rn_ishlat .AND. rn_ishlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  partial-slip' 
     165      ELSEIF ( 2. < rn_ishlat                      ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  strong-slip' 
    167166      ENDIF 
     167      ! 
     168      IF( .NOT. ln_landfast )   tau_icebfr(:,:) = 0._wp     ! NO Landfast ice : set to zero one for all 
    168169      ! 
    169170   END SUBROUTINE ice_rhg_init 
    170171 
     172#else 
     173   !!---------------------------------------------------------------------- 
     174   !!   Default option         Empty module          NO LIM-3 sea-ice model 
     175   !!---------------------------------------------------------------------- 
    171176#endif  
    172177 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90

    r8409 r8486  
    4545#  include "vectopt_loop_substitute.h90" 
    4646   !!---------------------------------------------------------------------- 
    47    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     47   !! NEMO/ICE 4.0 , NEMO Consortium (2017) 
    4848   !! $Id: icerhg_evp.F90 8378 2017-07-26 13:55:59Z clem $ 
    4949   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    386386         DO jj = 2, jpjm1 
    387387            DO ji = fs_2, fs_jpim1                
    388  
    389                ! U points 
     388               !                   !--- U points 
    390389               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             & 
    391390                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    & 
    … …  
    394393                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              & 
    395394                  &                  ) * r1_e1e2u(ji,jj) 
    396  
    397                ! V points 
     395                  ! 
     396                  !                !--- V points 
    398397               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             & 
    399398                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    & 
    … …  
    402401                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              & 
    403402                  &                  ) * r1_e1e2v(ji,jj) 
    404  
    405                ! u_ice at V point 
     403                  ! 
     404                  !                !--- u_ice at V point 
    406405               u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
    407406                  &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
    408                 
    409                ! v_ice at U point 
     407                  ! 
     408                  !                !--- v_ice at U point 
    410409               v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     & 
    411410                  &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
    412  
     411               ! 
    413412            END DO 
    414413         END DO 
    … …  
    417416         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp 
    418417         !  Bouillon et al. 2009 (eq 34-35) => stable 
    419          IF( MOD(jter,2) .EQ. 0 ) THEN ! even iterations 
    420              
     418         IF( MOD(jter,2) == 0 ) THEN ! even iterations 
     419            ! 
    421420            DO jj = 2, jpjm1 
    422421               DO ji = fs_2, fs_jpim1 
    423  
    424                   ! tau_io/(v_oce - v_ice) 
     422                     !                 !--- tau_io/(v_oce - v_ice) 
    425423                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
    426424                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 
    427  
    428                   ! Ocean-to-Ice stress 
     425                     !                 !--- Ocean-to-Ice stress 
    429426                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
    430  
    431                   ! tau_bottom/v_ice 
     427                     ! 
     428                     !                 !--- tau_bottom/v_ice 
    432429                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 
    433430                  zTauB = - tau_icebfr(ji,jj) / zvel 
    434  
    435                   ! Coriolis at V-points (energy conserving formulation) 
     431                     ! 
     432                     !                 !--- Coriolis at V-points (energy conserving formulation) 
    436433                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  & 
    437434                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
    438435                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
    439  
    440                   ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     436                     ! 
     437                     !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    441438                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
    442  
    443                   ! landfast switch => 0 = static friction ; 1 = sliding friction 
     439                     ! 
     440                     !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    444441                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
    445                    
    446                   ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
     442                     ! 
     443                     !                 !--- ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
    447444                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity 
    448445                     &                                     + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    … …  
    451448                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin 
    452449                     &           ) * zmaskV(ji,jj) 
     450                     ! 
    453451                  ! Bouillon 2013 
    454452                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  & 
    455453                  !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  & 
    456454                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj) 
    457  
     455                  ! 
    458456               END DO 
    459457            END DO 
    460458            CALL lbc_lnk( v_ice, 'V', -1. ) 
    461              
     459            ! 
    462460#if defined key_agrif 
    463461!!            CALL agrif_interp_lim3( 'V', jter, nn_nevp ) 
    … …  
    465463#endif 
    466464            IF( ln_bdy ) CALL bdy_ice_dyn( 'V' ) 
    467  
     465            ! 
    468466            DO jj = 2, jpjm1 
    469467               DO ji = fs_2, fs_jpim1 
    … …  
    505503            END DO 
    506504            CALL lbc_lnk( u_ice, 'U', -1. ) 
    507              
     505            ! 
    508506#if defined key_agrif 
    509507!!            CALL agrif_interp_lim3( 'U', jter, nn_nevp ) 
    … …  
    511509#endif 
    512510            IF( ln_bdy ) CALL bdy_ice_dyn( 'U' ) 
    513  
     511            ! 
    514512         ELSE ! odd iterations 
    515  
     513            ! 
    516514            DO jj = 2, jpjm1 
    517515               DO ji = fs_2, fs_jpim1 
    … …  
    553551            END DO 
    554552            CALL lbc_lnk( u_ice, 'U', -1. ) 
    555              
     553            ! 
    556554#if defined key_agrif 
    557555!!            CALL agrif_interp_lim3( 'U', jter, nn_nevp ) 
    … …  
    559557#endif 
    560558            IF( ln_bdy ) CALL bdy_ice_dyn( 'U' ) 
    561  
     559            ! 
    562560            DO jj = 2, jpjm1 
    563561               DO ji = fs_2, fs_jpim1 
    … …  
    599597            END DO 
    600598            CALL lbc_lnk( v_ice, 'V', -1. ) 
    601              
     599            ! 
    602600#if defined key_agrif 
    603601!!            CALL agrif_interp_lim3( 'V', jter, nn_nevp ) 
    … …  
    605603#endif 
    606604            IF( ln_bdy ) CALL bdy_ice_dyn( 'V' ) 
    607  
     605            ! 
    608606         ENDIF 
    609607          
    … …  
    675673      ! 5) SIMIP diagnostics 
    676674      !------------------------------------------------------------------------------! 
    677                             
     675 
     676!!gm  encapsulate with a flag (iom_use of the variable or better a flag defined one for all testing if one of the 
     677!!    diag is output.  NB the diag_...  are should only be ALLOCATED if the flag is true ! 
     678 
    678679      DO jj = 2, jpjm1 
    679680         DO ji = 2, jpim1 
    … …  
    714715 
    715716      CALL lbc_lnk_multi(   diag_sig1   , 'T',  1., diag_sig2   , 'T',  1.,   & 
    716                  &          diag_dssh_dx, 'U', -1., diag_dssh_dy, 'V', -1.,   & 
    717                  &          diag_corstrx, 'U', -1., diag_corstry, 'V', -1.,   &  
    718                  &          diag_intstrx, 'U', -1., diag_intstry, 'V', -1.    ) 
     717         &                  diag_dssh_dx, 'U', -1., diag_dssh_dy, 'V', -1.,   & 
     718         &                  diag_corstrx, 'U', -1., diag_corstry, 'V', -1.,   &  
     719         &                  diag_intstrx, 'U', -1., diag_intstry, 'V', -1.    ) 
    719720 
    720721      CALL lbc_lnk_multi(   diag_utau_oi, 'U', -1., diag_vtau_oi, 'V', -1.    ) 
    721722 
    722       CALL lbc_lnk_multi(   diag_xmtrp_ice, 'U', -1., diag_xmtrp_snw, 'U', -1., & 
    723                  &          diag_xatrp    , 'U', -1., diag_ymtrp_ice, 'V', -1., & 
    724                  &          diag_ymtrp_snw, 'V', -1., diag_yatrp    , 'V', -1.  ) 
     723      CALL lbc_lnk_multi(   diag_xmtrp_ice, 'U', -1., diag_xmtrp_snw, 'U', -1.,   & 
     724         &                  diag_xatrp    , 'U', -1., diag_ymtrp_ice, 'V', -1.,  & 
     725         &                  diag_ymtrp_snw, 'V', -1., diag_yatrp    , 'V', -1.    ) 
    725726 
    726727      ! 
    … …  
    734735         CALL prt_ctl_info(charout) 
    735736         CALL prt_ctl(tab2d_1=u_ice, clinfo1=' ice_rhg_evp  : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :') 
    736       ENDIF 
    737  
    738       ! print charge ellipse 
    739       ! This can be desactivated once the user is sure that the stress state 
    740       ! lie on the charge ellipse. See Bouillon et al. 08 for more details 
    741       IF(ln_ctl) THEN 
     737         ! 
     738         ! print charge ellipse 
     739         ! This can be desactivated once the user is sure that the stress state 
     740      ! lie on the charge ellipse. See Bouillon et al. (2008) for more details 
    742741         IF( MOD(kt_ice+nn_fsbc-1,nwrite) == 0 ) THEN 
    743742            WRITE(charout,FMT="('ice_rhg_evp  :', I4, I6, I1, I1, A10)") 1000, kt_ice, 0, 0, ' ch. ell. ' 
    … …  
    760759   END SUBROUTINE ice_rhg_evp 
    761760 
     761#else 
     762   !!---------------------------------------------------------------------- 
     763   !!   Default option         Empty module          NO LIM-3 sea-ice model 
     764   !!---------------------------------------------------------------------- 
    762765#endif 
    763766 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerst.F90

    r8413 r8486  
    44   !! Ice restart :  write the ice restart file 
    55   !!====================================================================== 
    6    !! History:   -   ! 2005-04 (M. Vancoppenolle) Original code 
    7    !!           3.0  ! 2008-03 (C. Ethe) restart files in using IOM interface 
    8    !!           4.0  ! 2011-02 (G. Madec) dynamical allocation 
     6   !! History:   3.0  ! 2005-04 (M. Vancoppenolle) Original code 
     7   !!             -   ! 2008-03 (C. Ethe) restart files in using IOM interface 
     8   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation 
    99   !!---------------------------------------------------------------------- 
    1010#if defined key_lim3 
    1111   !!---------------------------------------------------------------------- 
    12    !!   'key_lim3' :                                   LIM sea-ice model 
     12   !!   'key_lim3'                                        LIM sea-ice model 
    1313   !!---------------------------------------------------------------------- 
    1414   !!   ice_rst_opn   : open ice restart file 
    … …  
    3838 
    3939   !!---------------------------------------------------------------------- 
    40    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     40   !! NEMO/ICE 4.0 , NEMO Consortium (2017) 
    4141   !! $Id: icerst.F90 8411 2017-08-07 16:09:12Z clem $ 
    4242   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    123123      CALL iom_rstput( iter, nitrst, numriw, 'kt_ice' , REAL( iter   , wp ) )      ! date 
    124124 
     125!!gm   It is possible and easy to define a 3D domain size (jpi,jpj,jpl) or use a SIZE( tab, DIM=3) in iom_rtput ) 
     126!!gm         ===>>> just a simple   iom_rstput( iter, nitrst, numriw, 'v_i', v_i )  etc... 
     127!!gm   "just" ask Sebatien 
     128 
     129 
    125130      ! Prognostic variables  
    126131      DO jl = 1, jpl  
    … …  
    639644   !!   Default option :       Empty module            NO LIM sea-ice model 
    640645   !!---------------------------------------------------------------------- 
    641 CONTAINS 
    642    SUBROUTINE ice_rst_read             ! Empty routine 
    643    END SUBROUTINE ice_rst_read 
    644    SUBROUTINE ice_rst_write            ! Empty routine 
    645    END SUBROUTINE ice_rst_write 
    646646#endif 
    647647 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90

    r8426 r8486  
    1717#if defined key_lim3 
    1818   !!---------------------------------------------------------------------- 
    19    !!   'key_lim3' :                                  LIM 3.0 sea-ice model 
    20    !!---------------------------------------------------------------------- 
    21    !!   ice_stp  : sea-ice model time-stepping and update ocean sbc over ice-covered area 
    22    !!---------------------------------------------------------------------- 
    23    USE oce             ! ocean dynamics and tracers 
    24    USE dom_oce         ! ocean space and time domain 
    25    USE ice             ! LIM-3: ice variables 
    26    USE ice1D           ! LIM-3: thermodynamical variables 
     19   !!   'key_lim3'                                    LIM 3.0 sea-ice model 
     20   !!---------------------------------------------------------------------- 
     21   !!   ice_stp       : sea-ice model time-stepping and update ocean surf. boundary cond. over ice-covered area 
     22   !!   ice_init      : 
     23   !!   ice_run_init  :  
     24   !!---------------------------------------------------------------------- 
     25   USE oce            ! ocean dynamics and tracers 
     26   USE dom_oce        ! ocean space and time domain 
     27   USE c1d            ! 1D vertical configuration 
     28   USE ice            ! sea-ice: variables 
     29   USE ice1D          ! sea-ice: thermodynamical 1D variables 
    2730   ! 
    28    USE sbc_oce         ! Surface boundary condition: ocean fields 
    29    USE sbc_ice         ! Surface boundary condition: ice   fields 
    30    USE iceforcing      ! Surface boundary condition for sea ice 
     31   USE sbc_oce        ! Surface boundary condition: ocean fields 
     32   USE sbc_ice        ! Surface boundary condition: ice   fields 
     33   USE iceforcing     ! sea-ice: Surface boundary condition       !!gm why not icesbc module name 
    3134   ! 
    32    USE phycst          ! Define parameters for the routines 
    33    USE eosbn2          ! equation of state 
    34    USE icerhg          ! Ice rheology 
    35    USE iceadv          ! Ice advection 
    36    USE icethd          ! Ice thermodynamics 
    37    USE icerdgrft       ! Ice ridging/rafting 
    38    USE iceupdate       ! sea surface boundary condition 
    39    USE icedia          ! Ice budget diagnostics 
    40    USE icewri          ! Ice outputs 
    41    USE icerst          ! Ice restarts 
    42    USE icecor          ! Ice corrections 
    43    USE icevar          ! Ice variables switch 
    44    USE icectl          ! 
     35   USE phycst         ! Define parameters for the routines 
     36   USE eosbn2         ! equation of state 
     37   USE icerhg         ! sea-ice: rheology 
     38   USE iceadv         ! sea-ice: advection 
     39   USE icethd         ! sea-ice: thermodynamics 
     40   USE icerdgrft      ! sea-ice: ridging/rafting 
     41   USE iceupdate      ! sea-ice: sea surface boundary condition update 
     42   USE icedia         ! sea-ice: budget diagnostics 
     43   USE icewri         ! sea-ice: outputs 
     44   USE icerst         ! sea-ice: restarts 
     45   USE icecor         ! sea-ice: corrections 
     46   USE icevar         ! sea-ice: operations 
     47   USE icectl         ! sea-ice: control 
    4548   ! MV MP 2016 
    46    USE limmp 
     49   USE limmp          ! sea-ice: melt ponds 
    4750   ! END MV MP 2016 
    48    USE iceist          ! LIM initial state 
    49    USE icethd_sal      ! LIM ice thermodynamics: salinity 
     51   USE iceist         ! sea-ice: initial state 
     52   USE icethd_sal     ! sea-ice: thermodynamics and salinity 
    5053   ! 
    51    USE c1d             ! 1D vertical configuration 
    52    USE in_out_manager  ! I/O manager 
    53    USE iom             ! I/O manager library 
    54    USE prtctl          ! Print control 
    55    USE lib_fortran     ! 
    56    USE lbclnk          ! lateral boundary condition - MPP link 
    57    USE lib_mpp         ! MPP library 
    58    USE timing          ! Timing 
    59  
    60    USE bdy_oce   , ONLY: ln_bdy 
    61    USE bdyice          ! unstructured open boundary data 
     54   USE bdy_oce , ONLY : ln_bdy   ! flag for bdy 
     55   USE bdyice         ! unstructured open boundary data for sea-ice 
    6256# if defined key_agrif 
    6357   USE agrif_ice 
    … …  
    6559   USE agrif_lim3_interp 
    6660# endif 
     61   ! 
     62   USE in_out_manager ! I/O manager 
     63   USE iom            ! I/O manager library 
     64   USE prtctl         ! Print control 
     65   USE lib_fortran    !  
     66   USE lbclnk         ! lateral boundary condition - MPP link 
     67   USE lib_mpp        ! MPP library 
     68   USE timing         ! Timing 
    6769 
    6870   IMPLICIT NONE 
    6971   PRIVATE 
    7072 
    71    PUBLIC ice_stp  ! routine called by sbcmod.F90 
    72    PUBLIC ice_init ! routine called by sbcmod.F90 
     73   PUBLIC   ice_stp    ! called by sbcmod.F90 
     74   PUBLIC   ice_init   ! called by sbcmod.F90 
     75 
     76   INTEGER ::              nice_dyn   ! choice of the type of advection scheme 
     77   !                                  ! associated indices: 
     78   INTEGER, PARAMETER ::   np_dynNO   = 0   ! no ice dynamics and ice advection 
     79   INTEGER, PARAMETER ::   np_dynFULL = 1   ! full ice dynamics  (rheology + advection + ridging/rafting + correction) 
     80   INTEGER, PARAMETER ::   np_dyn     = 2   ! no ridging/rafting (rheology + advection                   + correction) 
     81   INTEGER, PARAMETER ::   np_dynPURE = 3   ! pure dynamics      (rheology + advection)  
    7382 
    7483   !! * Substitutions 
    7584#  include "vectopt_loop_substitute.h90" 
    7685   !!---------------------------------------------------------------------- 
    77    !! NEMO/OPA 4.0 , UCL NEMO Consortium (2011) 
     86   !! NEMO/ICE 4.0 , NEMO Consortium (2017) 
    7887   !! $Id: icestp.F90 8319 2017-07-11 15:00:44Z clem $ 
    7988   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    8594      !!                  ***  ROUTINE ice_stp  *** 
    8695      !! 
    87       !! ** Purpose :   update the ocean surface boundary condition via the 
    88       !!                Louvain la Neuve Sea Ice Model time stepping 
     96      !! ** Purpose :   sea-ice model time-stepping and update ocean surface 
     97      !!              boundary condition over ice-covered area 
    8998      !! 
    9099      !! ** Method  :   ice model time stepping 
    … …  
    102111      !!--------------------------------------------------------------------- 
    103112      INTEGER, INTENT(in) ::   kt      ! ocean time step 
    104       INTEGER, INTENT(in) ::   ksbc    ! type of sbc flux ( 1 = user defined formulation,  
    105                                        !                    3 = bulk formulation, 
    106                                        !                    4 = Pure Coupled formulation) 
    107       INTEGER  ::   jl                 ! dummy loop index 
    108       !!---------------------------------------------------------------------- 
    109  
     113      INTEGER, INTENT(in) ::   ksbc    ! flux formulation (user defined, bulk, or Pure Coupled) 
     114      ! 
     115      INTEGER ::   jl   ! dummy loop index 
     116      !!---------------------------------------------------------------------- 
     117      ! 
    110118      IF( nn_timing == 1 )   CALL timing_start('ice_stp') 
    111  
    112       !-----------------------! 
    113       ! --- Ice time step --- ! 
    114       !-----------------------! 
    115       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 
    116           
    117          ! Ice model time step 
    118          kt_ice = kt 
    119  
     119      ! 
     120      !                                      !-----------------------! 
     121      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN   ! --- Ice time step --- ! 
     122         !                                   !-----------------------! 
     123         ! 
     124         kt_ice = kt       ! Ice model time step 
     125         ! 
    120126# if defined key_agrif 
    121127         IF( .NOT. Agrif_Root() )  lim_nbstep = MOD( lim_nbstep, Agrif_irhot() * Agrif_Parent(nn_fsbc) / nn_fsbc ) + 1 
    122128# endif 
    123129 
    124          ! mean surface ocean current at ice velocity point (C-grid dynamics :  U- & V-points as the ocean) 
     130         !                 ! mean surface ocean current at ice velocity point 
    125131         u_oce(:,:) = ssu_m(:,:) * umask(:,:,1) 
    126132         v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1) 
    127  
    128          ! masked sea surface freezing temperature [Kelvin] (set to rt0 over land) 
     133!!gm the umask, vmask above is useless as ssu_m, ssv_m are build from masked un,vn (used t be here for B-grid sea-ice) 
     134 
     135         !                 ! masked sea surface freezing temperature [Kelvin] (set to rt0 over land) 
    129136         CALL eos_fzp( sss_m(:,:) , t_bo(:,:) ) 
    130137         t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) ) 
    131138         ! 
    132                                       CALL ice_bef         ! Store previous ice values 
     139         CALL ice_bef      ! Store previous ice values 
     140         ! 
    133141         !------------------------------------------------! 
    134142         ! --- Dynamical coupling with the atmosphere --- ! 
    135143         !------------------------------------------------! 
    136          ! it provides: 
     144         ! It provides the following fields used in sea ice model: 
    137145         !    utau_ice, vtau_ice = surface ice stress [N/m2] 
    138          !-------------------------------------------------- 
    139                                       CALL ice_forcing_tau( kt, ksbc, utau_ice, vtau_ice ) 
     146         !------------------------------------------------! 
     147         CALL ice_forcing_tau( kt, ksbc, utau_ice, vtau_ice ) 
    140148                                       
    141          !-------------------------------------------------------! 
    142          ! --- ice dynamics and transport (except in 1D case) ---! 
    143          !-------------------------------------------------------! 
    144                                       CALL ice_diag0           ! set diag of mass, heat and salt fluxes to 0 
    145                                       CALL ice_rst_opn( kt )   ! Open Ice restart file 
    146          ! 
    147          ! --- zap this if no ice dynamics --- ! 
    148          IF( .NOT. lk_c1d .AND. ln_limdyn ) THEN 
    149                                       CALL ice_rhg( kt )       ! -- rheology   
    150                                       CALL ice_adv( kt )       ! -- advection 
    151             IF( nn_limdyn == 2 .AND. nn_monocat /= 2 )  &      ! -- ridging/rafting 
    152                &                      CALL ice_rdgrft( kt )          
    153             IF( nn_limdyn == 2 )      CALL ice_cor( kt , 1 )   ! -- Corrections 
    154             ! 
    155          ENDIF 
    156          ! --- 
    157           
     149         !-------------------------------------! 
     150         ! --- ice dynamics and advection  --- ! 
     151         !-------------------------------------! 
     152         CALL ice_diag0             ! set diag of mass, heat and salt fluxes to 0 
     153         CALL ice_rst_opn( kt )     ! Open Ice restart file (if necessary)  
     154         ! 
     155         SELECT CASE( nice_dyn )    
     156         CASE ( np_dynFULL )          !==  all dynamical processes  ==! 
     157            CALL ice_rhg   ( kt )          ! -- rheology   
     158            CALL ice_adv   ( kt )          ! -- advection of ice 
     159            CALL ice_rdgrft( kt )          ! -- ridging/rafting  
     160            CALL ice_cor   ( kt , 1 )      ! -- Corrections 
     161         CASE ( np_dyn )              !==  pure dynamics only ==!   (no ridging/rafting)   (nono cat. case 2) 
     162            CALL ice_rhg   ( kt )          ! -- rheology   
     163            CALL ice_adv   ( kt )          ! -- advection of ice 
     164            CALL ice_cor   ( kt , 1 )      ! -- Corrections 
     165         CASE ( np_dynPURE )          !==  pure dynamics only ==!   (nn_limdyn= 0 or 1 ) 
     166            CALL ice_rhg   ( kt )          ! -- rheology   
     167            CALL ice_adv   ( kt )          ! -- advection of ice 
     168         END SELECT 
     169         ! 
     170         !                          !==  lateral boundary conditions  ==! 
    158171#if defined key_agrif 
    159          IF( .NOT. Agrif_Root() )     CALL agrif_interp_lim3('T') 
     172         IF( .NOT. Agrif_Root()     )   CALL agrif_interp_lim3('T')   ! -- AGRIF  
    160173#endif 
    161          IF( ln_limthd .AND. ln_bdy ) CALL bdy_ice( kt )       ! -- bdy ice thermo  
    162          ! previous lead fraction and ice volume for flux calculations 
    163                                       CALL ice_var_glo2eqv     ! ht_i and ht_s for ice albedo calculation 
    164                                       CALL ice_var_agg(1)      ! at_i for coupling  
    165                                       CALL ice_bef 
     174         IF( ln_limthd .AND. ln_bdy )   CALL bdy_ice( kt )            ! -- bdy ice thermo 
     175         ! 
     176         ! 
     177         !                          !==  previous lead fraction and ice volume for flux calculations 
     178         ! 
     179         CALL ice_var_glo2eqv            ! ht_i and ht_s for ice albedo calculation 
     180         CALL ice_var_agg(1)             ! at_i for coupling  
     181         CALL ice_bef                    ! Store previous ice values 
    166182 
    167183         !------------------------------------------------------! 
    … …  
    169185         !------------------------------------------------------! 
    170186         ! It provides the following fields used in sea ice model: 
    171          !    fr1_i0  , fr2_i0                         = 1sr & 2nd fraction of qsr penetration in ice  [%] 
    172          !    emp_oce , emp_ice                        = E-P over ocean and sea ice                    [Kg/m2/s] 
    173          !    sprecip                                  = solid precipitation                           [Kg/m2/s] 
    174          !    evap_ice                                 = sublimation                                   [Kg/m2/s] 
    175          !    qsr_tot , qns_tot                        = solar & non solar heat flux (total)           [W/m2] 
    176          !    qsr_ice , qns_ice                        = solar & non solar heat flux over ice          [W/m2] 
    177          !    dqns_ice                                 = non solar  heat sensistivity                  [W/m2] 
    178          !    qemp_oce, qemp_ice, qprec_ice, qevap_ice = sensible heat (associated with evap & precip) [W/m2] 
    179          !------------------------------------------------------------------------------------------------------ 
    180                                       CALL ice_forcing_flx( kt, ksbc ) 
     187         !    fr1_i0  , fr2_i0     = 1sr & 2nd fraction of qsr penetration in ice  [%] 
     188         !    emp_oce , emp_ice    = E-P over ocean and sea ice                    [Kg/m2/s] 
     189         !    sprecip              = solid precipitation                           [Kg/m2/s] 
     190         !    evap_ice             = sublimation                                   [Kg/m2/s] 
     191         !    qsr_tot , qns_tot    = solar & non solar heat flux (total)           [W/m2] 
     192         !    qsr_ice , qns_ice    = solar & non solar heat flux over ice          [W/m2] 
     193         !    dqns_ice             = non solar  heat sensistivity                  [W/m2] 
     194         !    qemp_oce, qemp_ice,  = sensible heat (associated with evap & precip) [W/m2] 
     195         !    qprec_ice, qevap_ice 
     196         !------------------------------------------------------! 
     197         CALL ice_forcing_flx( kt, ksbc ) 
    181198 
    182199         !----------------------------! 
    183200         ! --- ice thermodynamics --- ! 
    184201         !----------------------------! 
    185          ! --- zap this if no ice thermo --- ! 
    186          IF( ln_limthd )              CALL ice_thd( kt )        ! -- Ice thermodynamics       
     202         IF( ln_limthd )            CALL ice_thd( kt )          ! -- Ice thermodynamics       
    187203 
    188204         ! MV MP 2016 
    189          IF ( ln_pnd )                CALL lim_mp( kt )         ! -- Melt ponds 
     205         IF ( ln_pnd )              CALL lim_mp( kt )           ! -- Melt ponds 
    190206         ! END MV MP 2016 
    191207 
    192          IF( ln_limthd )              CALL ice_cor( kt , 2 )    ! -- Corrections 
     208         IF( ln_limthd )            CALL ice_cor( kt , 2 )      ! -- Corrections 
    193209         ! --- 
    194210# if defined key_agrif 
    195          IF( .NOT. Agrif_Root() )     CALL agrif_update_lim3( kt ) 
     211         IF( .NOT. Agrif_Root() )   CALL agrif_update_lim3( kt ) 
    196212# endif 
    197                                       CALL ice_var_glo2eqv      ! necessary calls (at least for coupling) 
    198                                       CALL ice_var_agg( 2 )     ! necessary calls (at least for coupling) 
    199                                       ! 
     213                                    CALL ice_var_glo2eqv        ! necessary calls (at least for coupling) 
     214                                    CALL ice_var_agg( 2 )       ! necessary calls (at least for coupling) 
     215                                    ! 
    200216!! clem: one should switch the calculation of the fluxes onto the parent grid but the following calls do not work 
    201217!!       moreover it should only be called at the update frequency (as in agrif_lim3_update.F90) 
    … …  
    203219!!         IF( .NOT. Agrif_Root() )   CALL Agrif_ChildGrid_To_ParentGrid() 
    204220!!# endif 
    205                                       CALL ice_update_flx( kt )    ! -- Update surface ocean mass, heat and salt fluxes 
     221                                    CALL ice_update_flx( kt )   ! -- Update ocean surface mass, heat and salt fluxes 
    206222!!# if defined key_agrif 
    207223!!         IF( .NOT. Agrif_Root() )   CALL Agrif_ParentGrid_To_ChildGrid() 
    208224!!# endif 
    209          IF( ln_limdiahsb )           CALL ice_dia( kt )     ! -- Diagnostics and outputs  
     225         IF( ln_limdiahsb )           CALL ice_dia( kt )        ! -- Diagnostics and outputs  
    210226         ! 
    211227                                      CALL ice_wri( 1 )         ! -- Ice outputs  
    212228         ! 
    213          IF( kt == nit000 .AND. ln_rstart )   & 
     229         IF( kt == nit000 .AND. ln_rstart )   &                !!gm  I don't understand the ln_rstart, if needed, add a comment, please 
    214230            &                         CALL iom_close( numrir )  ! close input ice restart file 
    215231         ! 
    … …  
    287303      IF( ln_limdiahsb) CALL ice_dia_init  ! initialization for diags 
    288304      ! 
    289       fr_i(:,:)     = at_i(:,:)         ! initialisation of sea-ice fraction 
     305      fr_i  (:,:)   = at_i(:,:)         ! initialisation of sea-ice fraction 
    290306      tn_ice(:,:,:) = t_su(:,:,:)       ! initialisation of surface temp for coupled simu 
    291307      ! 
    292308      DO jj = 1, jpj 
    293309         DO ji = 1, jpi 
    294             IF( gphit(ji,jj) > 0._wp ) THEN  ;  rn_amax_2d(ji,jj) = rn_amax_n  ! NH 
    295             ELSE                             ;  rn_amax_2d(ji,jj) = rn_amax_s  ! SH 
     310            IF( gphit(ji,jj) > 0._wp ) THEN   ;   rn_amax_2d(ji,jj) = rn_amax_n  ! NH 
     311            ELSE                              ;   rn_amax_2d(ji,jj) = rn_amax_s  ! SH 
    296312            ENDIF 
    297313         END DO 
    … …  
    318334      !!------------------------------------------------------------------- 
    319335      ! 
    320       REWIND( numnam_ice_ref )              ! Namelist namicerun in reference namelist : Parameters for ice 
     336      REWIND( numnam_ice_ref )      ! Namelist namicerun in reference namelist : Parameters for ice 
    321337      READ  ( numnam_ice_ref, namicerun, IOSTAT = ios, ERR = 901) 
    322338901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in reference namelist', lwp ) 
    323339 
    324       REWIND( numnam_ice_cfg )              ! Namelist namicerun in configuration namelist : Parameters for ice 
     340      REWIND( numnam_ice_cfg )      ! Namelist namicerun in configuration namelist : Parameters for ice 
    325341      READ  ( numnam_ice_cfg, namicerun, IOSTAT = ios, ERR = 902 ) 
    326342902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in configuration namelist', lwp ) 
    327343      IF(lwm) WRITE ( numoni, namicerun ) 
    328344      ! 
    329       REWIND( numnam_ice_ref )              ! Namelist namicediag in reference namelist : Parameters for ice 
     345      REWIND( numnam_ice_ref )      ! Namelist namicediag in reference namelist : Parameters for ice 
    330346      READ  ( numnam_ice_ref, namicediag, IOSTAT = ios, ERR = 903) 
    331347903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicediag in reference namelist', lwp ) 
    332348 
    333       REWIND( numnam_ice_cfg )              ! Namelist namicediag in configuration namelist : Parameters for ice 
     349      REWIND( numnam_ice_cfg )      ! Namelist namicediag in configuration namelist : Parameters for ice 
    334350      READ  ( numnam_ice_cfg, namicediag, IOSTAT = ios, ERR = 904 ) 
    335351904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicediag in configuration namelist', lwp ) 
    336352      IF(lwm) WRITE ( numoni, namicediag ) 
    337353      ! 
    338       IF(lwp) THEN                        ! control print 
     354      IF(lwp) THEN                  ! control print 
    339355         WRITE(numout,*) 
    340356         WRITE(numout,*) 'ice_run_init : ice share parameters for dynamics/advection/thermo of sea-ice' 
    341357         WRITE(numout,*) ' ~~~~~~' 
    342          WRITE(numout,*) '   number of ice  categories                              jpl    = ', jpl 
    343          WRITE(numout,*) '   number of ice  layers                                  nlay_i = ', nlay_i 
    344          WRITE(numout,*) '   number of snow layers                                  nlay_s = ', nlay_s 
    345          WRITE(numout,*) '   virtual ITD mono-category param (1-4) or not (0)   nn_monocat = ', nn_monocat 
    346          WRITE(numout,*) '   maximum ice concentration for NH                              = ', rn_amax_n  
    347          WRITE(numout,*) '   maximum ice concentration for SH                              = ', rn_amax_s 
    348          WRITE(numout,*) '   Ice thermodynamics (T) or not (F)                  ln_limthd  = ', ln_limthd 
    349          WRITE(numout,*) '   Ice dynamics       (T) or not (F)                  ln_limdyn  = ', ln_limdyn 
    350          WRITE(numout,*) '     (ln_limdyn=T) Ice dynamics switch                nn_limdyn  = ', nn_limdyn 
    351          WRITE(numout,*) '       2: total' 
    352          WRITE(numout,*) '       1: advection only (no diffusion, no ridging/rafting)' 
    353          WRITE(numout,*) '       0: advection only (as 1 + prescribed velocity, bypass rheology)' 
    354          WRITE(numout,*) '     (ln_limdyn=T) prescribed u-vel (case nn_limdyn=0)           = ', rn_uice 
    355          WRITE(numout,*) '     (ln_limdyn=T) prescribed v-vel (case nn_limdyn=0)           = ', rn_vice 
     358         WRITE(numout,*) '   Namelist namicerun : ' 
     359         WRITE(numout,*) '      number of ice  categories                              jpl    = ', jpl 
     360         WRITE(numout,*) '      number of ice  layers                                  nlay_i = ', nlay_i 
     361         WRITE(numout,*) '      number of snow layers                                  nlay_s = ', nlay_s 
     362         WRITE(numout,*) '      virtual ITD mono-category param (1-4) or not (0)   nn_monocat = ', nn_monocat 
     363         WRITE(numout,*) '      maximum ice concentration for NH                              = ', rn_amax_n  
     364         WRITE(numout,*) '      maximum ice concentration for SH                              = ', rn_amax_s 
     365         WRITE(numout,*) '      Ice thermodynamics (T) or not (F)                  ln_limthd  = ', ln_limthd 
     366         WRITE(numout,*) '      Ice dynamics       (T) or not (F)                  ln_limdyn  = ', ln_limdyn 
     367         WRITE(numout,*) '         associated switch                               nn_limdyn  = ', nn_limdyn 
     368         WRITE(numout,*) '            =2 all processes (default option)' 
     369         WRITE(numout,*) '            =1 advection only (no ridging/rafting)' 
     370         WRITE(numout,*) '            =0 advection only with prescribed velocity given by ' 
     371         WRITE(numout,*) '               a uniform field        (u,v)_ice = (rn_uice,rn_vice) = (', rn_uice,',', rn_vice,')' 
    356372         WRITE(numout,*) 
    357          WRITE(numout,*) '...and ice diagnostics' 
    358          WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~' 
    359          WRITE(numout,*) '   Diagnose online heat/mass/salt budget     ln_limdiachk  = ', ln_limdiachk 
    360          WRITE(numout,*) '   Output          heat/mass/salt budget     ln_limdiahsb  = ', ln_limdiahsb 
    361          WRITE(numout,*) '   control prints in ocean.out for (i,j)=(iiceprt,jiceprt) = ', ln_limctl 
    362          WRITE(numout,*) '   i-index for control prints (ln_limctl=true)             = ', iiceprt 
    363          WRITE(numout,*) '   j-index for control prints (ln_limctl=true)             = ', jiceprt 
     373         WRITE(numout,*) '   Namelist namicediag : ' 
     374         WRITE(numout,*) '      Diagnose online heat/mass/salt budget      ln_limdiachk = ', ln_limdiachk 
     375         WRITE(numout,*) '      Output          heat/mass/salt budget      ln_limdiahsb = ', ln_limdiahsb 
     376         WRITE(numout,*) '      control prints for a given grid point         ln_limctl = ', ln_limctl 
     377         WRITE(numout,*) '         chosen grid point position         (iiceprt,jiceprt) = (', iiceprt,',', jiceprt,')' 
    364378      ENDIF 
    365379      ! 
    366       IF ( ( jpl > 1 ) .AND. ( nn_monocat == 1 ) ) THEN 
     380      !                                         !--- check consistency 
     381      IF ( jpl > 1 .AND. nn_monocat == 1 ) THEN 
    367382         nn_monocat = 0 
    368383         IF(lwp) WRITE(numout,*) 
    369384         IF(lwp) WRITE(numout,*) '   nn_monocat forced to 0 as jpl>1, i.e. multi-category case is chosen' 
    370385      ENDIF 
    371       IF ( ( jpl == 1 ) .AND. ( nn_monocat == 0 ) ) THEN 
     386      IF ( jpl == 1 .AND. nn_monocat == 0 ) THEN 
    372387         CALL ctl_stop( 'STOP', 'ice_run_init : if jpl=1 then nn_monocat should be between 1 and 4' ) 
    373388      ENDIF 
    374389      ! 
    375       ! sea-ice timestep and inverse 
    376       rdt_ice   = REAL(nn_fsbc) * rdt   
    377       r1_rdtice = 1._wp / rdt_ice  
    378  
    379       ! inverse of nlay_i and nlay_s 
    380       r1_nlay_i = 1._wp / REAL( nlay_i, wp ) 
     390      IF( ln_bdy .AND. ln_limdiachk )   CALL ctl_warn('ice_run_init: online conservation check does not work with BDY') 
     391      ! 
     392      !                             ! set the choice of ice dynamics 
     393      IF( lk_c1d .OR. .NOT.ln_limdyn ) THEN 
     394         nice_dyn = np_dynNO                    !--- no dynamics 
     395      ELSE 
     396         SELECT CASE( nn_limdyn ) 
     397         CASE( 2 )                     
     398            IF( nn_monocat /= 2 ) THEN          !--- full dynamics (rheology + advection + ridging/rafting + correction) 
     399               nice_dyn = np_dynFULL 
     400            ELSE 
     401               nice_dyn = np_dyn                !--- dynamics without ridging/rafting 
     402            ENDIF 
     403         CASE( 0 , 1 )                          !--- dynamics without ridging/rafting and correction  
     404            nice_dyn = np_dynPURE 
     405         END SELECT 
     406      ENDIF 
     407      !                                         !--- simple conservative piling, comparable with LIM2 
     408      l_piling = nn_limdyn == 1 .OR. ( nn_monocat == 2  .AND.  jpl == 1 ) 
     409      ! 
     410      rdt_ice   = REAL(nn_fsbc) * rdt           !--- sea-ice timestep and inverse 
     411      r1_rdtice = 1._wp / rdt_ice 
     412      IF( lwp ) WRITE(numout,*) '   ice timestep rdt_ice  = ', rdt_ice 
     413      ! 
     414      r1_nlay_i = 1._wp / REAL( nlay_i, wp )    !--- inverse of nlay_i and nlay_s 
    381415      r1_nlay_s = 1._wp / REAL( nlay_s, wp ) 
    382       ! 
    383       IF( lwp .AND. ln_bdy .AND. ln_limdiachk )  & 
    384          &   CALL ctl_warn('online conservation check activated but it does not work with BDY') 
    385       ! 
    386       IF( lwp ) WRITE(numout,*) '   ice timestep rdt_ice  = ', rdt_ice 
    387416      ! 
    388417   END SUBROUTINE ice_run_init 
    … …  
    397426      !! ** input   :   Namelist namiceitd 
    398427      !!------------------------------------------------------------------- 
    399       INTEGER  ::   ios                 ! Local integer output status for namelist read 
     428      INTEGER  ::   jl    ! dummy loop index 
     429      INTEGER  ::   ios   ! Local integer output status for namelist read 
     430      REAL(wp) ::   zc1, zc2, zc3, zx1          ! local scalars 
     431      REAL(wp) ::   zhmax, znum, zden, zalpha   !   -      - 
     432      !! 
    400433      NAMELIST/namiceitd/ rn_himean 
    401       ! 
    402       INTEGER  ::   jl                   ! dummy loop index 
    403       REAL(wp) ::   zc1, zc2, zc3, zx1   ! local scalars 
    404       REAL(wp) ::   zhmax, znum, zden, zalpha ! 
    405434      !!------------------------------------------------------------------ 
    406435      ! 
    407       REWIND( numnam_ice_ref )              ! Namelist namiceitd in reference namelist : Parameters for ice 
     436      REWIND( numnam_ice_ref )      ! Namelist namiceitd in reference namelist : Parameters for ice 
    408437      READ  ( numnam_ice_ref, namiceitd, IOSTAT = ios, ERR = 905) 
    409438905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in reference namelist', lwp ) 
    410439 
    411       REWIND( numnam_ice_cfg )              ! Namelist namiceitd in configuration namelist : Parameters for ice 
     440      REWIND( numnam_ice_cfg )      ! Namelist namiceitd in configuration namelist : Parameters for ice 
    412441      READ  ( numnam_ice_cfg, namiceitd, IOSTAT = ios, ERR = 906 ) 
    413442906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in configuration namelist', lwp ) 
    414443      IF(lwm) WRITE ( numoni, namiceitd ) 
    415444      ! 
    416       IF(lwp) THEN                        ! control print 
     445      IF(lwp) THEN                  ! control print 
    417446         WRITE(numout,*) 
    418447         WRITE(numout,*) 'ice_itd_init : Initialization of ice cat distribution ' 
    419448         WRITE(numout,*) '~~~~~~~~~~~~' 
    420          WRITE(numout,*) '   mean ice thickness in the domain            rn_himean = ', rn_himean 
     449         WRITE(numout,*) '   Namelist namicerun : ' 
     450         WRITE(numout,*) '      mean ice thickness in the domain               rn_himean = ', rn_himean 
    421451      ENDIF 
    422452      ! 
    423       !---------------------------------- 
    424       !- Thickness categories boundaries 
    425       !---------------------------------- 
    426       ! 
    427       !==  h^(-alpha) function  ==! 
    428       zalpha = 0.05_wp 
     453      !-----------------------------------! 
     454      !  Thickness categories boundaries  ! 
     455      !-----------------------------------! 
     456      ! 
     457      zalpha = 0.05_wp              ! max of each category (from h^(-alpha) function) 
    429458      zhmax  = 3._wp * rn_himean 
    430459      DO jl = 1, jpl 
    … …  
    441470      ! 
    442471      IF(lwp) WRITE(numout,*) 
    443       IF(lwp) WRITE(numout,*) '      Thickness category boundaries ' 
    444       IF(lwp) WRITE(numout,*) '         hi_max ', hi_max(0:jpl) 
     472      IF(lwp) WRITE(numout,*) '   ===>>>   resulting thickness category boundaries :' 
     473      IF(lwp) WRITE(numout,*) '            hi_max(:)= ', hi_max(0:jpl) 
    445474      ! 
    446475   END SUBROUTINE ice_itd_init 
     476 
    447477 
    448478   SUBROUTINE ice_bef 
    … …  
    472502         END DO    
    473503