Changeset 8515


Ignore:
Timestamp:
2017-09-08T18:19:17+02:00 (3 years ago)
Author:
clem
Message:

changes in style - part5 - very nearly finished

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

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref

    r8514 r8515  
    11!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    2 !! LIM3 namelist:   
     2!! ESIM namelist:   
    33!!              1 - Generic parameters                 (namice_run) 
    44!!              2 - Ice thickness discretization       (namice_itd) 
     
    3535!------------------------------------------------------------------------------ 
    3636   rn_himean        =   2.0           !  expected domain-average ice thickness (m) 
     37   rn_himin         =   0.1           !  minimum ice thickness (m) used in remapping 
    3738/ 
    3839!------------------------------------------------------------------------------ 
     
    4647      rn_perdg      =   17.0          !     ridging work divided by pot. energy change in ridging 
    4748                   ! -- ice_rdgrft -- ! 
    48    rn_cs            =   0.5           !  fraction of shearing energy contributing to ridging 
     49   rn_csrdg         =   0.5           !  fraction of shearing energy contributing to ridging 
    4950              ! -- ice_rdgrft_prep -- ! 
    5051   ln_partf_lin     = .false.         !  Linear ridging participation function (Thorndike et al, 1975) 
     
    100101!------------------------------------------------------------------------------ 
    101102   ln_icethd        = .true.            !  ice thermo   (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 
    102                      ! -- limthd_dif -- ! 
     103                     ! -- icethd_dif -- ! 
    103104   rn_kappa_i       =   1.0             !  radiation attenuation coefficient in sea ice (m-1) 
    104105   ln_cndi_U64      = .false.           !  sea ice thermal conductivity: k = k0 + beta.S/T            (Untersteiner, 1964) 
     
    107108   rn_cnd_s         =   0.31            !  thermal conductivity of the snow (0.31 W/m/K, Maykut and Untersteiner, 1971) 
    108109                                        !     Obs: 0.1-0.5 (Lecomte et al, JAMES 2013) 
    109                       ! -- limthd_dh -- ! 
     110                      ! -- icethd_dh -- ! 
    110111   ln_icedH         = .true.            !  activate ice thickness change from growing/melting (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 
    111112   rn_blow_s        =   0.66            !     mesure of snow blowing into the leads 
    112113                                        !     = 1 => no snow blowing, < 1 => some snow blowing 
    113                       ! -- limthd_da -- ! 
     114                      ! -- icethd_da -- ! 
    114115   ln_icedA         = .true.            !  activate lateral melting param. (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 
    115116      rn_beta       =   1.0             !     coef. beta for lateral melting param. Recommended range=[0.8-1.2] 
     
    120121                                        !      => 6  vs 8m = +40% melting at the peak (A~0.5) 
    121122                                        !         10 vs 8m = -20% melting 
    122                      ! -- limthd_lac -- ! 
     123                     ! -- icethd_lac -- ! 
    123124   ln_icedO         = .true.            !  activate ice growth in open-water (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 
    124       rn_hinew      =   0.1             !     thickness for new ice formation in open water (m) 
     125      rn_hinew      =   0.1             !     thickness for new ice formation in open water (m), must be larger than rn_hnewice 
    125126      ln_frazil     = .false.           !     Frazil ice parameterization (ice collection as a function of wind) 
    126127         rn_maxfraz =   1.0             !        maximum fraction of frazil ice collecting at the ice base 
    127128         rn_vfraz   =   0.417           !        thresold drift speed for frazil ice collecting at the ice bottom (m/s) 
    128129         rn_Cfraz   =   5.0             !        squeezing coefficient for frazil ice collecting at the ice bottom 
    129                     ! -- limitd_th -- ! 
    130    rn_himin         =   0.1             !  minimum ice thickness (m) used in remapping, must be smaller than rn_hnewice 
    131130                         ! -- icestp -- ! 
    132131   nn_iceflx        =  -1               !  Redistribute heat flux over ice categories 
     
    143142&namice_sal     !   Ice salinity 
    144143!------------------------------------------------------------------------------ 
    145                     ! -- limthd_sal -- ! 
     144                    ! -- icethd_sal -- ! 
    146145   ln_icedS       = .true.             !  activate gravity drainage and flushing (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 
    147146   nn_icesal      =   2                !  ice salinity option 
     
    171170&namice_ini     !   Ice initialization 
    172171!------------------------------------------------------------------------------ 
    173                      ! -- limistate -- ! 
     172                     ! -- iceistate -- ! 
    174173   ln_iceini      = .true.             !  activate ice initialization (T) or not (F) 
    175174   ln_iceini_file = .false.            !  netcdf file provided for initialization (T) or not (F) 
    176    rn_thres_sst   =   2.0              !  maximum water temperature with initial ice (degC) 
     175   rn_thres_sst   =   2.0              !  max delta temp. above Tfreeze with initial ice = (sst - tfreeze) 
    177176   rn_hts_ini_n   =   0.3              !  initial real snow thickness (m), North 
    178177   rn_hts_ini_s   =   0.3              !        "            "             South 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90

    r8514 r8515  
    3636   PRIVATE 
    3737 
    38    PUBLIC   ice_istate      ! called by icestp.F90 
     38   PUBLIC   ice_istate        ! called by icestp.F90 
     39   PUBLIC   ice_istate_init   ! called by icestp.F90 
    3940 
    4041   INTEGER , PARAMETER ::   jpfldi = 6           ! maximum number of files to read 
     
    102103 
    103104      IF(lwp) WRITE(numout,*) 
    104       IF(lwp) WRITE(numout,*) 'ice_istate : sea-ice initialization ' 
    105       IF(lwp) WRITE(numout,*) '~~~~~~~~~~ ' 
     105      IF(lwp) WRITE(numout,*) 'ice_istate: sea-ice initialization ' 
     106      IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
    106107 
    107108      !-------------------------------------------------------------------- 
     
    109110      !-------------------------------------------------------------------- 
    110111      ! 
    111       CALL ice_istate_init 
    112  
    113112      ! init surface temperature 
    114113      DO jl = 1, jpl 
     
    154153            !----------------------------- 
    155154            ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array 
    156             DO jj = 1, jpj 
    157                DO ji = 1, jpi 
    158                   IF( ff_t(ji,jj) >= 0._wp ) THEN 
    159                      zht_i_ini(ji,jj) = rn_hti_ini_n * zswitch(ji,jj) 
    160                      zht_s_ini(ji,jj) = rn_hts_ini_n * zswitch(ji,jj) 
    161                      zat_i_ini(ji,jj) = rn_ati_ini_n * zswitch(ji,jj) 
    162                      zts_u_ini(ji,jj) = rn_tmi_ini_n * zswitch(ji,jj) 
    163                      zsm_i_ini(ji,jj) = rn_smi_ini_n * zswitch(ji,jj) 
    164                      ztm_i_ini(ji,jj) = rn_tmi_ini_n * zswitch(ji,jj) 
    165                   ELSE 
    166                      zht_i_ini(ji,jj) = rn_hti_ini_s * zswitch(ji,jj) 
    167                      zht_s_ini(ji,jj) = rn_hts_ini_s * zswitch(ji,jj) 
    168                      zat_i_ini(ji,jj) = rn_ati_ini_s * zswitch(ji,jj) 
    169                      zts_u_ini(ji,jj) = rn_tmi_ini_s * zswitch(ji,jj) 
    170                      zsm_i_ini(ji,jj) = rn_smi_ini_s * zswitch(ji,jj) 
    171                      ztm_i_ini(ji,jj) = rn_tmi_ini_s * zswitch(ji,jj) 
    172                   ENDIF 
    173                END DO 
    174             END DO 
     155            WHERE( ff_t(:,:) >= 0._wp ) 
     156               zht_i_ini(:,:) = rn_hti_ini_n * zswitch(:,:) 
     157               zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:) 
     158               zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:) 
     159               zts_u_ini(:,:) = rn_tmi_ini_n * zswitch(:,:) 
     160               zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:) 
     161               ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:) 
     162            ELSEWHERE 
     163               zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:) 
     164               zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:) 
     165               zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:) 
     166               zts_u_ini(:,:) = rn_tmi_ini_s * zswitch(:,:) 
     167               zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:) 
     168               ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:) 
     169            END WHERE 
    175170            ! 
    176171         ENDIF ! ln_iceini_file 
     
    554549      slf_i(jp_ati) = sn_ati  ;  slf_i(jp_tsu) = sn_tsu 
    555550      slf_i(jp_tmi) = sn_tmi  ;  slf_i(jp_smi) = sn_smi 
    556  
    557       ! Define the initial parameters 
    558       ! ------------------------- 
    559  
    560       IF(lwp) THEN 
     551      ! 
     552      ! 
     553      IF(lwp) THEN                          ! control print 
    561554         WRITE(numout,*) 
    562555         WRITE(numout,*) 'ice_istate_init: ice parameters inititialisation ' 
    563556         WRITE(numout,*) '~~~~~~~~~~~~~~~' 
    564557         WRITE(numout,*) '   Namelist namice_ini' 
    565          WRITE(numout,*) '      initialization with ice (T) or not (F)       ln_iceini     = ', ln_iceini 
    566          WRITE(numout,*) '      ice initialization from a netcdf file      ln_iceini_file  = ', ln_iceini_file 
    567          WRITE(numout,*) '      threshold water temp. for initial sea-ice    rn_thres_sst  = ', rn_thres_sst 
    568          WRITE(numout,*) '      initial snow thickness in the north          rn_hts_ini_n  = ', rn_hts_ini_n 
    569          WRITE(numout,*) '      initial snow thickness in the south          rn_hts_ini_s  = ', rn_hts_ini_s  
    570          WRITE(numout,*) '      initial ice thickness  in the north          rn_hti_ini_n  = ', rn_hti_ini_n 
    571          WRITE(numout,*) '      initial ice thickness  in the south          rn_hti_ini_s  = ', rn_hti_ini_s 
    572          WRITE(numout,*) '      initial ice concentr.  in the north          rn_ati_ini_n  = ', rn_ati_ini_n 
    573          WRITE(numout,*) '      initial ice concentr.  in the north          rn_ati_ini_s  = ', rn_ati_ini_s 
    574          WRITE(numout,*) '      initial  ice salinity  in the north          rn_smi_ini_n  = ', rn_smi_ini_n 
    575          WRITE(numout,*) '      initial  ice salinity  in the south          rn_smi_ini_s  = ', rn_smi_ini_s 
    576          WRITE(numout,*) '      initial  ice/snw temp  in the north          rn_tmi_ini_n  = ', rn_tmi_ini_n 
    577          WRITE(numout,*) '      initial  ice/snw temp  in the south          rn_tmi_ini_s  = ', rn_tmi_ini_s 
     558         WRITE(numout,*) '      initialization with ice (T) or not (F)                 ln_iceini     = ', ln_iceini 
     559         WRITE(numout,*) '      ice initialization from a netcdf file                ln_iceini_file  = ', ln_iceini_file 
     560         WRITE(numout,*) '      max delta ocean temp. above Tfreeze with initial ice   rn_thres_sst  = ', rn_thres_sst 
     561         WRITE(numout,*) '      initial snow thickness in the north                    rn_hts_ini_n  = ', rn_hts_ini_n 
     562         WRITE(numout,*) '      initial snow thickness in the south                    rn_hts_ini_s  = ', rn_hts_ini_s  
     563         WRITE(numout,*) '      initial ice thickness  in the north                    rn_hti_ini_n  = ', rn_hti_ini_n 
     564         WRITE(numout,*) '      initial ice thickness  in the south                    rn_hti_ini_s  = ', rn_hti_ini_s 
     565         WRITE(numout,*) '      initial ice concentr.  in the north                    rn_ati_ini_n  = ', rn_ati_ini_n 
     566         WRITE(numout,*) '      initial ice concentr.  in the north                    rn_ati_ini_s  = ', rn_ati_ini_s 
     567         WRITE(numout,*) '      initial  ice salinity  in the north                    rn_smi_ini_n  = ', rn_smi_ini_n 
     568         WRITE(numout,*) '      initial  ice salinity  in the south                    rn_smi_ini_s  = ', rn_smi_ini_s 
     569         WRITE(numout,*) '      initial  ice/snw temp  in the north                    rn_tmi_ini_n  = ', rn_tmi_ini_n 
     570         WRITE(numout,*) '      initial  ice/snw temp  in the south                    rn_tmi_ini_s  = ', rn_tmi_ini_s 
    578571      ENDIF 
    579572 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90

    r8514 r8515  
    650650      REAL(wp) ::   zhmax, znum, zden, zalpha   !   -      - 
    651651      !! 
    652       NAMELIST/namice_itd/ rn_himean 
     652      NAMELIST/namice_itd/ rn_himean, rn_himin 
    653653      !!------------------------------------------------------------------ 
    654654      ! 
     
    668668         WRITE(numout,*) '   Namelist namice_itd : ' 
    669669         WRITE(numout,*) '      mean ice thickness in the domain               rn_himean = ', rn_himean 
     670         WRITE(numout,*) '      minimum ice thickness                          rn_himin  = ', rn_himin  
    670671      ENDIF 
    671672      ! 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90

    r8514 r8515  
    5252   ! 
    5353   ! ** namelist (namice_rdgrft) ** 
    54    REAL(wp) ::   rn_cs            ! fraction of shearing energy contributing to ridging             
     54   REAL(wp) ::   rn_csrdg         ! fraction of shearing energy contributing to ridging             
    5555   LOGICAL  ::   ln_partf_lin     ! participation function linear (Thorndike et al. (1975)) 
    5656   REAL(wp) ::   rn_gstar         !    fractional area of young ice contributing to ridging 
     
    170170            !  (thick, newly ridged ice). 
    171171 
    172             closing_net(ji,jj) = rn_cs * 0.5_wp * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp ) 
     172            closing_net(ji,jj) = rn_csrdg * 0.5_wp * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp ) 
    173173 
    174174            ! 2.2 divu_adv 
     
    867867      NAMELIST/namice_rdgrft/ ln_str_H79, rn_pstar, rn_crhg, & 
    868868         &                    ln_str_R75, rn_perdg,          & 
    869          &                    rn_cs     ,                    & 
     869         &                    rn_csrdg  ,                    & 
    870870         &                    ln_partf_lin, rn_gstar,        & 
    871871         &                    ln_partf_exp, rn_astar,        &  
     
    893893         WRITE(numout,*) '      ice strength parameterization Rothrock (1975)            ln_str_R75   = ', ln_str_R75  
    894894         WRITE(numout,*) '            Ratio of ridging work to PotEner change in ridging rn_perdg     = ', rn_perdg  
    895          WRITE(numout,*) '      Fraction of shear energy contributing to ridging         rn_cs        = ', rn_cs  
     895         WRITE(numout,*) '      Fraction of shear energy contributing to ridging         rn_csrdg     = ', rn_csrdg  
    896896         WRITE(numout,*) '      linear ridging participation function                    ln_partf_lin = ', ln_partf_lin 
    897897         WRITE(numout,*) '            Fraction of ice coverage contributing to ridging   rn_gstar     = ', rn_gstar 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90

    r8512 r8515  
    5252CONTAINS 
    5353 
    54    SUBROUTINE ice_rhg_evp( kt, stress1_i, stress2_i, stress12_i, u_ice, v_ice, shear_i, divu_i, delta_i ) 
     54   SUBROUTINE ice_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, pu_ice, pv_ice, pshear_i, pdivu_i, pdelta_i ) 
    5555      !!------------------------------------------------------------------- 
    5656      !!                 ***  SUBROUTINE ice_rhg_evp  *** 
     
    108108      !!------------------------------------------------------------------- 
    109109      INTEGER, INTENT(in) ::   kt      ! time step 
    110       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   stress1_i, stress2_i, stress12_i 
    111       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   u_ice, v_ice, shear_i, divu_i, delta_i  
     110      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i 
     111      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pu_ice, pv_ice, pshear_i, pdivu_i, pdelta_i  
    112112      !! 
    113113      INTEGER ::   ji, jj       ! dummy loop indices 
     
    247247 
    248248      ! Initialise stress tensor  
    249       zs1 (:,:) = stress1_i (:,:)  
    250       zs2 (:,:) = stress2_i (:,:) 
    251       zs12(:,:) = stress12_i(:,:) 
     249      zs1 (:,:) = pstress1_i (:,:)  
     250      zs2 (:,:) = pstress2_i (:,:) 
     251      zs12(:,:) = pstress12_i(:,:) 
    252252 
    253253      ! Ice strength 
     
    340340         IF(ln_ctl) THEN   ! Convergence test 
    341341            DO jj = 1, jpjm1 
    342                zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 
    343                zv_ice(:,jj) = v_ice(:,jj) 
     342               zu_ice(:,jj) = pu_ice(:,jj) ! velocity at previous time step 
     343               zv_ice(:,jj) = pv_ice(:,jj) 
    344344            END DO 
    345345         ENDIF 
     
    350350 
    351351               ! shear at F points 
    352                zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   & 
    353                   &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
     352               zds(ji,jj) = ( ( pu_ice(ji,jj+1) * r1_e1u(ji,jj+1) - pu_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   & 
     353                  &         + ( pv_ice(ji+1,jj) * r1_e2v(ji+1,jj) - pv_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
    354354                  &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) 
    355355 
     
    367367               
    368368               ! divergence at T points 
    369                zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
    370                   &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
     369               zdiv  = ( e2u(ji,jj) * pu_ice(ji,jj) - e2u(ji-1,jj) * pu_ice(ji-1,jj)   & 
     370                  &    + e1v(ji,jj) * pv_ice(ji,jj) - e1v(ji,jj-1) * pv_ice(ji,jj-1)   & 
    371371                  &    ) * r1_e1e2t(ji,jj) 
    372372               zdiv2 = zdiv * zdiv 
    373373                
    374374               ! tension at T points 
    375                zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
    376                   &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
     375               zdt  = ( ( pu_ice(ji,jj) * r1_e2u(ji,jj) - pu_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
     376                  &   - ( pv_ice(ji,jj) * r1_e1v(ji,jj) - pv_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
    377377                  &   ) * r1_e1e2t(ji,jj) 
    378378               zdt2 = zdt * zdt 
     
    426426                  ! 
    427427                  !                !--- u_ice at V point 
    428                u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
    429                   &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
     428               u_iceV(ji,jj) = 0.5_wp * ( ( pu_ice(ji,jj  ) + pu_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
     429                  &                     + ( pu_ice(ji,jj+1) + pu_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
    430430                  ! 
    431431                  !                !--- v_ice at U point 
    432                v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     & 
    433                   &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
     432               v_iceU(ji,jj) = 0.5_wp * ( ( pv_ice(ji  ,jj) + pv_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     & 
     433                  &                     + ( pv_ice(ji+1,jj) + pv_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
    434434               ! 
    435435            END DO 
     
    444444               DO ji = fs_2, fs_jpim1 
    445445                     !                 !--- tau_io/(v_oce - v_ice) 
    446                   zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
     446                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( pv_ice(ji,jj) - v_oce (ji,jj) ) * ( pv_ice(ji,jj) - v_oce (ji,jj) )  & 
    447447                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 
    448448                     !                 !--- Ocean-to-Ice stress 
    449                   ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
     449                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - pv_ice(ji,jj) ) 
    450450                     ! 
    451451                     !                 !--- tau_bottom/v_ice 
    452                   zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 
     452                  zvel  = MAX( zepsi, SQRT( pv_ice(ji,jj) * pv_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 
    453453                  zTauB = - tau_icebfr(ji,jj) / zvel 
    454454                     ! 
    455455                     !                 !--- Coriolis at V-points (energy conserving formulation) 
    456456                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  & 
    457                      &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
    458                      &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
     457                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * pu_ice(ji,jj  ) + e2u(ji-1,jj  ) * pu_ice(ji-1,jj  ) )  & 
     458                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * pu_ice(ji,jj+1) + e2u(ji-1,jj+1) * pu_ice(ji-1,jj+1) ) ) 
    459459                     ! 
    460460                     !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     
    465465                     ! 
    466466                     !                 !--- ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
    467                   v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity 
    468                      &                                     + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    469                      &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast 
    470                      &             + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
    471                      &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin 
    472                      &           ) * zmaskV(ji,jj) 
     467                  pv_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * pv_ice(ji,jj)                            &  ! previous velocity 
     468                     &                                      + zTauE + zTauO * pv_ice(ji,jj)                            &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     469                     &                                      ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast 
     470                     &              + ( 1._wp - rswitch ) * pv_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
     471                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                   &  ! v_ice = v_oce if mass < zmmin 
     472                     &            ) * zmaskV(ji,jj) 
    473473                     ! 
    474474                  ! Bouillon 2013 
    475                   !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  & 
     475                  !!pv_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * pv_ice(ji,jj) + pv_ice_b(ji,jj) )                  & 
    476476                  !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  & 
    477477                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj) 
     
    479479               END DO 
    480480            END DO 
    481             CALL lbc_lnk( v_ice, 'V', -1. ) 
     481            CALL lbc_lnk( pv_ice, 'V', -1. ) 
    482482            ! 
    483483#if defined key_agrif 
     
    491491                                
    492492                  ! tau_io/(u_oce - u_ice) 
    493                   zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
     493                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( pu_ice(ji,jj) - u_oce (ji,jj) ) * ( pu_ice(ji,jj) - u_oce (ji,jj) )  & 
    494494                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
    495495 
    496496                  ! Ocean-to-Ice stress 
    497                   ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
     497                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - pu_ice(ji,jj) ) 
    498498 
    499499                  ! tau_bottom/u_ice 
    500                   zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) ) 
     500                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + pu_ice(ji,jj) * pu_ice(ji,jj) ) ) 
    501501                  zTauB = - tau_icebfr(ji,jj) / zvel 
    502502 
    503503                  ! Coriolis at U-points (energy conserving formulation) 
    504504                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  & 
    505                      &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
    506                      &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
     505                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * pv_ice(ji  ,jj) + e1v(ji  ,jj-1) * pv_ice(ji  ,jj-1) )  & 
     506                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * pv_ice(ji+1,jj) + e1v(ji+1,jj-1) * pv_ice(ji+1,jj-1) ) ) 
    507507                   
    508508                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     
    513513 
    514514                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
    515                   u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity 
    516                      &                                     + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    517                      &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast 
    518                      &             + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
    519                      &             ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin  
    520                      &           ) * zmaskU(ji,jj) 
     515                  pu_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * pu_ice(ji,jj)                            &  ! previous velocity 
     516                     &                                      + zTauE + zTauO * pu_ice(ji,jj)                            &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     517                     &                                      ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast 
     518                     &              + ( 1._wp - rswitch ) * pu_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
     519                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                   &  ! v_ice = v_oce if mass < zmmin  
     520                     &            ) * zmaskU(ji,jj) 
    521521                  ! Bouillon 2013 
    522                   !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  & 
     522                  !!pu_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * pu_ice(ji,jj) + pu_ice_b(ji,jj) )                  & 
    523523                  !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  & 
    524524                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj) 
    525525               END DO 
    526526            END DO 
    527             CALL lbc_lnk( u_ice, 'U', -1. ) 
     527            CALL lbc_lnk( pu_ice, 'U', -1. ) 
    528528            ! 
    529529#if defined key_agrif 
     
    539539 
    540540                  ! tau_io/(u_oce - u_ice) 
    541                   zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
     541                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( pu_ice(ji,jj) - u_oce (ji,jj) ) * ( pu_ice(ji,jj) - u_oce (ji,jj) )  & 
    542542                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
    543543 
    544544                  ! Ocean-to-Ice stress 
    545                   ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
     545                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - pu_ice(ji,jj) ) 
    546546 
    547547                  ! tau_bottom/u_ice 
    548                   zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) ) 
     548                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + pu_ice(ji,jj) * pu_ice(ji,jj) ) ) 
    549549                  zTauB = - tau_icebfr(ji,jj) / zvel 
    550550 
    551551                  ! Coriolis at U-points (energy conserving formulation) 
    552552                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  & 
    553                      &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
    554                      &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
     553                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * pv_ice(ji  ,jj) + e1v(ji  ,jj-1) * pv_ice(ji  ,jj-1) )  & 
     554                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * pv_ice(ji+1,jj) + e1v(ji+1,jj-1) * pv_ice(ji+1,jj-1) ) ) 
    555555 
    556556                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     
    561561 
    562562                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
    563                   u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity 
    564                      &                                     + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    565                      &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast 
    566                      &             + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
    567                      &             ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin 
    568                      &           ) * zmaskU(ji,jj) 
     563                  pu_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * pu_ice(ji,jj)                            &  ! previous velocity 
     564                     &                                      + zTauE + zTauO * pu_ice(ji,jj)                            &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     565                     &                                      ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast 
     566                     &              + ( 1._wp - rswitch ) * pu_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
     567                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                   &  ! v_ice = v_oce if mass < zmmin 
     568                     &            ) * zmaskU(ji,jj) 
    569569                  ! Bouillon 2013 
    570                   !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  & 
     570                  !!pu_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * pu_ice(ji,jj) + pu_ice_b(ji,jj) )                  & 
    571571                  !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  & 
    572572                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj) 
    573573               END DO 
    574574            END DO 
    575             CALL lbc_lnk( u_ice, 'U', -1. ) 
     575            CALL lbc_lnk( pu_ice, 'U', -1. ) 
    576576            ! 
    577577#if defined key_agrif 
     
    585585          
    586586                  ! tau_io/(v_oce - v_ice) 
    587                   zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
     587                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( pv_ice(ji,jj) - v_oce (ji,jj) ) * ( pv_ice(ji,jj) - v_oce (ji,jj) )  & 
    588588                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 
    589589 
    590590                  ! Ocean-to-Ice stress 
    591                   ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
     591                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - pv_ice(ji,jj) ) 
    592592 
    593593                  ! tau_bottom/v_ice 
    594                   zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 
     594                  zvel  = MAX( zepsi, SQRT( pv_ice(ji,jj) * pv_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 
    595595                  ztauB = - tau_icebfr(ji,jj) / zvel 
    596596                   
    597597                  ! Coriolis at V-points (energy conserving formulation) 
    598598                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  & 
    599                      &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
    600                      &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
     599                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * pu_ice(ji,jj  ) + e2u(ji-1,jj  ) * pu_ice(ji-1,jj  ) )  & 
     600                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * pu_ice(ji,jj+1) + e2u(ji-1,jj+1) * pu_ice(ji-1,jj+1) ) ) 
    601601 
    602602                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     
    607607                   
    608608                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
    609                   v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity 
    610                      &                                     + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    611                      &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast 
    612                      &             + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
    613                      &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin 
    614                      &           ) * zmaskV(ji,jj) 
     609                  pv_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * pv_ice(ji,jj)                            &  ! previous velocity 
     610                     &                                      + zTauE + zTauO * pv_ice(ji,jj)                            &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     611                     &                                      ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast 
     612                     &              + ( 1._wp - rswitch ) * pv_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
     613                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                   &  ! v_ice = v_oce if mass < zmmin 
     614                     &            ) * zmaskV(ji,jj) 
    615615                  ! Bouillon 2013 
    616                   !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  & 
     616                  !!pv_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * pv_ice(ji,jj) + pv_ice_b(ji,jj) )                  & 
    617617                  !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  & 
    618618                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj) 
    619619               END DO 
    620620            END DO 
    621             CALL lbc_lnk( v_ice, 'V', -1. ) 
     621            CALL lbc_lnk( pv_ice, 'V', -1. ) 
    622622            ! 
    623623#if defined key_agrif 
     
    631631         IF(ln_ctl) THEN   ! Convergence test 
    632632            DO jj = 2 , jpjm1 
    633                zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 
     633               zresr(:,jj) = MAX( ABS( pu_ice(:,jj) - zu_ice(:,jj) ), ABS( pv_ice(:,jj) - zv_ice(:,jj) ) ) 
    634634            END DO 
    635635            zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) ) 
     
    648648 
    649649            ! shear at F points 
    650             zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   & 
    651                &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
     650            zds(ji,jj) = ( ( pu_ice(ji,jj+1) * r1_e1u(ji,jj+1) - pu_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   & 
     651               &         + ( pv_ice(ji+1,jj) * r1_e2v(ji+1,jj) - pv_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
    652652               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) 
    653653 
     
    660660             
    661661            ! tension**2 at T points 
    662             zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
    663                &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
     662            zdt  = ( ( pu_ice(ji,jj) * r1_e2u(ji,jj) - pu_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
     663               &   - ( pv_ice(ji,jj) * r1_e1v(ji,jj) - pv_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
    664664               &   ) * r1_e1e2t(ji,jj) 
    665665            zdt2 = zdt * zdt 
     
    671671             
    672672            ! shear at T points 
    673             shear_i(ji,jj) = SQRT( zdt2 + zds2 ) 
     673            pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) 
    674674 
    675675            ! divergence at T points 
    676             divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
    677                &            + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
    678                &            ) * r1_e1e2t(ji,jj) 
     676            pdivu_i(ji,jj) = ( e2u(ji,jj) * pu_ice(ji,jj) - e2u(ji-1,jj) * pu_ice(ji-1,jj)   & 
     677               &             + e1v(ji,jj) * pv_ice(ji,jj) - e1v(ji,jj-1) * pv_ice(ji,jj-1)   & 
     678               &             ) * r1_e1e2t(ji,jj) 
    679679             
    680680            ! delta at T points 
    681             zdelta         = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )   
     681            zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )   
    682682            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0 
    683             delta_i(ji,jj) = zdelta + rn_creepl * rswitch 
     683            pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch 
    684684 
    685685         END DO 
    686686      END DO 
    687       CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1., delta_i, 'T', 1. ) 
     687      CALL lbc_lnk_multi( pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. ) 
    688688       
    689689      ! --- Store the stress tensor for the next time step --- ! 
    690       stress1_i (:,:) = zs1 (:,:) 
    691       stress2_i (:,:) = zs2 (:,:) 
    692       stress12_i(:,:) = zs12(:,:) 
     690      pstress1_i (:,:) = zs1 (:,:) 
     691      pstress2_i (:,:) = zs2 (:,:) 
     692      pstress12_i(:,:) = zs12(:,:) 
    693693      ! 
    694694 
     
    703703 
    704704      ! --- divergence, shear and strength --- ! 
    705       IF( iom_use('idive')  )   CALL iom_put( "idive"  , divu_i(:,:)   * zswi(:,:) )   ! divergence 
    706       IF( iom_use('ishear') )   CALL iom_put( "ishear" , shear_i(:,:)  * zswi(:,:) )   ! shear 
     705      IF( iom_use('idive')  )   CALL iom_put( "idive"  , pdivu_i(:,:)   * zswi(:,:) )   ! divergence 
     706      IF( iom_use('ishear') )   CALL iom_put( "ishear" , pshear_i(:,:)  * zswi(:,:) )   ! shear 
    707707      IF( iom_use('icestr') )   CALL iom_put( "icestr" , strength(:,:) * zswi(:,:) )   ! Ice strength 
    708708 
     
    714714         DO jj = 2, jpjm1 
    715715            DO ji = 2, jpim1 
    716                zdum1 = ( zswi(ji-1,jj) * stress12_i(ji-1,jj) + zswi(ji  ,jj-1) * stress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point 
    717                   &     zswi(ji  ,jj) * stress12_i(ji  ,jj) + zswi(ji-1,jj-1) * stress12_i(ji-1,jj-1) )  & 
    718                   &   / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) ) 
    719  
    720                zshear = SQRT( stress2_i(ji,jj) * stress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress   
     716               zdum1 = ( zswi(ji-1,jj) * pstress12_i(ji-1,jj) + zswi(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point 
     717                  &      zswi(ji  ,jj) * pstress12_i(ji  ,jj) + zswi(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  & 
     718                  &    / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) ) 
     719 
     720               zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress   
    721721 
    722722               zdum2 = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 
    723723 
    724 !!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 
    725 !!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) 
    726 !!               zsig3(ji,jj) = zdum2**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress 
     724!!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 
     725!!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) 
     726!!               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress 
    727727!!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 
    728                zsig1(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015 
     728               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015 
    729729               zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                    ! shear stress 
    730                zsig3(ji,jj) = zdum2**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 
     730               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 
    731731            END DO 
    732732         END DO 
     
    773773                
    774774               ! 2D ice mass, snow mass, area transport arrays (X, Y) 
    775                zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch 
    776                zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch 
     775               zfac_x = 0.5 * pu_ice(ji,jj) * e2u(ji,jj) * rswitch 
     776               zfac_y = 0.5 * pv_ice(ji,jj) * e1v(ji,jj) * rswitch 
    777777                
    778778               zdiag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90

    r8514 r8515  
    285285      !                                ! Initial sea-ice state 
    286286      IF( .NOT. ln_rstart ) THEN              ! start from rest: sea-ice deduced from sst 
     287         CALL ice_istate_init 
    287288         CALL ice_istate 
    288289      ELSE                                    ! start from a restart file 
     
    301302      tn_ice(:,:,:) = t_su(:,:,:)       ! initialisation of surface temp for coupled simu 
    302303      ! 
    303       DO jj = 1, jpj 
    304          DO ji = 1, jpi 
    305             IF( gphit(ji,jj) > 0._wp ) THEN   ;   rn_amax_2d(ji,jj) = rn_amax_n  ! NH 
    306             ELSE                              ;   rn_amax_2d(ji,jj) = rn_amax_s  ! SH 
    307             ENDIF 
    308          END DO 
    309       END DO 
     304      WHERE( gphit(:,:) > 0._wp )   ;   rn_amax_2d(:,:) = rn_amax_n  ! NH 
     305      ELSEWHERE                     ;   rn_amax_2d(:,:) = rn_amax_s  ! SH 
     306      END WHERE 
    310307      ! 
    311308   END SUBROUTINE ice_init 
     
    381378      !!---------------------------------------------------------------------- 
    382379      ! 
    383       DO jl = 1, jpl 
    384          DO jj = 2, jpjm1 
    385             DO ji = 2, jpim1 
    386                a_i_b  (ji,jj,jl)   = a_i  (ji,jj,jl)     ! ice area 
    387                v_i_b  (ji,jj,jl)   = v_i  (ji,jj,jl)     ! ice volume 
    388                v_s_b  (ji,jj,jl)   = v_s  (ji,jj,jl)     ! snow volume 
    389                smv_i_b(ji,jj,jl)   = smv_i(ji,jj,jl)     ! salt content 
    390                oa_i_b (ji,jj,jl)   = oa_i (ji,jj,jl)     ! areal age content 
    391                e_s_b  (ji,jj,:,jl) = e_s  (ji,jj,:,jl)   ! snow thermal energy 
    392                e_i_b  (ji,jj,:,jl) = e_i  (ji,jj,:,jl)   ! ice thermal energy 
    393                !                                         ! ice thickness 
    394                rswitch = MAX( 0._wp , SIGN( 1._wp, a_i_b(ji,jj,jl) - epsi20 ) ) ! 0 if no ice and 1 if yes 
    395                ht_i_b(ji,jj,jl) = v_i_b (ji,jj,jl) / MAX( a_i_b(ji,jj,jl) , epsi20 ) * rswitch 
    396                ht_s_b(ji,jj,jl) = v_s_b (ji,jj,jl) / MAX( a_i_b(ji,jj,jl) , epsi20 ) * rswitch 
    397             END DO 
    398          END DO    
    399       END DO 
    400       CALL lbc_lnk_multi(  a_i_b, 'T', 1., v_i_b , 'T', 1., ht_i_b, 'T', 1., smv_i_b, 'T', 1.,   & 
    401          &                oa_i_b, 'T', 1., v_s_b , 'T', 1., ht_s_b, 'T', 1. ) 
    402       CALL lbc_lnk( e_i_b, 'T', 1. ) 
    403       CALL lbc_lnk( e_s_b, 'T', 1. ) 
     380      a_i_b  (:,:,:)   = a_i  (:,:,:)     ! ice area 
     381      v_i_b  (:,:,:)   = v_i  (:,:,:)     ! ice volume 
     382      v_s_b  (:,:,:)   = v_s  (:,:,:)     ! snow volume 
     383      smv_i_b(:,:,:)   = smv_i(:,:,:)     ! salt content 
     384      oa_i_b (:,:,:)   = oa_i (:,:,:)     ! areal age content 
     385      e_s_b  (:,:,:,:) = e_s  (:,:,:,:)   ! snow thermal energy 
     386      e_i_b  (:,:,:,:) = e_i  (:,:,:,:)   ! ice thermal energy 
     387      WHERE( a_i_b(:,:,:) >= epsi20 ) 
     388         ht_i_b(:,:,:) = v_i_b (:,:,:) / a_i_b(:,:,:)   ! ice thickness 
     389         ht_s_b(:,:,:) = v_s_b (:,:,:) / a_i_b(:,:,:)   ! snw thickness 
     390      ELSEWHERE 
     391         ht_i_b(:,:,:) = 0._wp 
     392         ht_s_b(:,:,:) = 0._wp 
     393      END WHERE 
    404394       
    405 !!gm Question:  here , a_i_b, u_ice and v_ice  are defined over the whole domain,  
    406 !!gm            so why not just a copy over the whole domain and no lbc_lnk ? 
    407 !!gm            that is some thing like: 
    408 !            at_i_b(:,:)  = SUM( a_i_b(:,:,:), dim=3 ) 
    409 !            u_ice_b(:,:) = u_ice(:,:) 
    410 !            v_ice_b(:,:) = v_ice(:,:) 
    411 !    idem for the loop above 
    412 !!gm 
    413395      ! ice velocities & total concentration 
    414       DO jj = 2, jpjm1 
    415          DO ji = 2, jpim1 
    416             at_i_b(ji,jj)  = SUM( a_i_b(ji,jj,:) ) 
    417             u_ice_b(ji,jj) = u_ice(ji,jj) 
    418             v_ice_b(ji,jj) = v_ice(ji,jj) 
    419          END DO 
    420       END DO 
    421       CALL lbc_lnk_multi( at_i_b, 'T', 1., u_ice_b , 'U', -1., v_ice_b , 'V', -1. ) 
     396      at_i_b(:,:)  = SUM( a_i_b(:,:,:), dim=3 ) 
     397      u_ice_b(:,:) = u_ice(:,:) 
     398      v_ice_b(:,:) = v_ice(:,:) 
    422399      ! 
    423400   END SUBROUTINE store_fields 
     
    433410      INTEGER  ::   ji, jj      ! dummy loop index 
    434411      !!---------------------------------------------------------------------- 
    435       DO jj = 1, jpj 
    436          DO ji = 1, jpi 
    437             sfx    (ji,jj) = 0._wp   ; 
    438             sfx_bri(ji,jj) = 0._wp   ;   sfx_lam(ji,jj) = 0._wp 
    439             sfx_sni(ji,jj) = 0._wp   ;   sfx_opw(ji,jj) = 0._wp 
    440             sfx_bog(ji,jj) = 0._wp   ;   sfx_dyn(ji,jj) = 0._wp 
    441             sfx_bom(ji,jj) = 0._wp   ;   sfx_sum(ji,jj) = 0._wp 
    442             sfx_res(ji,jj) = 0._wp   ;   sfx_sub(ji,jj) = 0._wp 
    443             ! 
    444             wfx_snw(ji,jj) = 0._wp   ;   wfx_ice(ji,jj) = 0._wp 
    445             wfx_sni(ji,jj) = 0._wp   ;   wfx_opw(ji,jj) = 0._wp 
    446             wfx_bog(ji,jj) = 0._wp   ;   wfx_dyn(ji,jj) = 0._wp 
    447             wfx_bom(ji,jj) = 0._wp   ;   wfx_sum(ji,jj) = 0._wp 
    448             wfx_res(ji,jj) = 0._wp   ;   wfx_sub(ji,jj) = 0._wp 
    449             wfx_spr(ji,jj) = 0._wp   ;   wfx_lam(ji,jj) = 0._wp   
    450             wfx_snw_dyn(ji,jj) = 0._wp ; wfx_snw_sum(ji,jj) = 0._wp 
    451             wfx_snw_sub(ji,jj) = 0._wp ; wfx_ice_sub(ji,jj) = 0._wp 
    452             wfx_snw_sni(ji,jj) = 0._wp  
    453             ! MV MP 2016 
    454             wfx_pnd(ji,jj) = 0._wp 
    455             ! END MV MP 2016 
    456              
    457             hfx_thd(ji,jj) = 0._wp   ; 
    458             hfx_snw(ji,jj) = 0._wp   ;   hfx_opw(ji,jj) = 0._wp 
    459             hfx_bog(ji,jj) = 0._wp   ;   hfx_dyn(ji,jj) = 0._wp 
    460             hfx_bom(ji,jj) = 0._wp   ;   hfx_sum(ji,jj) = 0._wp 
    461             hfx_res(ji,jj) = 0._wp   ;   hfx_sub(ji,jj) = 0._wp 
    462             hfx_spr(ji,jj) = 0._wp   ;   hfx_dif(ji,jj) = 0._wp 
    463             hfx_err(ji,jj) = 0._wp   ;   hfx_err_rem(ji,jj) = 0._wp 
    464             hfx_err_dif(ji,jj) = 0._wp 
    465             wfx_err_sub(ji,jj) = 0._wp 
    466             ! 
    467             afx_tot(ji,jj) = 0._wp   ; 
    468             ! 
    469             diag_heat(ji,jj) = 0._wp ;   diag_smvi(ji,jj) = 0._wp 
    470             diag_vice(ji,jj) = 0._wp ;   diag_vsnw(ji,jj) = 0._wp 
    471              
    472             ! SIMIP diagnostics 
    473             diag_fc_bo(ji,jj)    = 0._wp ; diag_fc_su(ji,jj)    = 0._wp 
    474              
    475             tau_icebfr(ji,jj) = 0._wp; ! landfast ice param only (clem: important to keep the init here) 
    476          END DO 
    477       END DO 
     412      sfx    (:,:) = 0._wp   ; 
     413      sfx_bri(:,:) = 0._wp   ;   sfx_lam(:,:) = 0._wp 
     414      sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp 
     415      sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp 
     416      sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp 
     417      sfx_res(:,:) = 0._wp   ;   sfx_sub(:,:) = 0._wp 
     418      ! 
     419      wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp 
     420      wfx_sni(:,:) = 0._wp   ;   wfx_opw(:,:) = 0._wp 
     421      wfx_bog(:,:) = 0._wp   ;   wfx_dyn(:,:) = 0._wp 
     422      wfx_bom(:,:) = 0._wp   ;   wfx_sum(:,:) = 0._wp 
     423      wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp 
     424      wfx_spr(:,:) = 0._wp   ;   wfx_lam(:,:) = 0._wp   
     425      wfx_snw_dyn(:,:) = 0._wp ; wfx_snw_sum(:,:) = 0._wp 
     426      wfx_snw_sub(:,:) = 0._wp ; wfx_ice_sub(:,:) = 0._wp 
     427      wfx_snw_sni(:,:) = 0._wp  
     428      ! MV MP 2016 
     429      wfx_pnd(:,:) = 0._wp 
     430      ! END MV MP 2016 
     431 
     432      hfx_thd(:,:) = 0._wp   ; 
     433      hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp 
     434      hfx_bog(:,:) = 0._wp   ;   hfx_dyn(:,:) = 0._wp 
     435      hfx_bom(:,:) = 0._wp   ;   hfx_sum(:,:) = 0._wp 
     436      hfx_res(:,:) = 0._wp   ;   hfx_sub(:,:) = 0._wp 
     437      hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp 
     438      hfx_err(:,:) = 0._wp   ;   hfx_err_rem(:,:) = 0._wp 
     439      hfx_err_dif(:,:) = 0._wp 
     440      wfx_err_sub(:,:) = 0._wp 
     441      ! 
     442      afx_tot(:,:) = 0._wp   ; 
     443      ! 
     444      diag_heat(:,:) = 0._wp ;   diag_smvi(:,:) = 0._wp 
     445      diag_vice(:,:) = 0._wp ;   diag_vsnw(:,:) = 0._wp 
     446 
     447      ! SIMIP diagnostics 
     448      diag_fc_bo(:,:)    = 0._wp ; diag_fc_su(:,:)    = 0._wp 
     449 
     450      tau_icebfr(:,:) = 0._wp; ! landfast ice param only (clem: important to keep the init here) 
    478451       
    479452   END SUBROUTINE ice_diag0 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90

    r8514 r8515  
    539539      INTEGER  ::   ios   ! Local integer output status for namelist read 
    540540      !! 
    541       NAMELIST/namice_thd/ ln_icethd, rn_kappa_i, ln_cndi_U64, ln_cndi_P07, ln_dqns_i, rn_cnd_s,          & 
    542          &                ln_icedH, rn_blow_s,                                                             & 
    543          &                ln_icedA, rn_beta, rn_dmin,                                                     & 
    544          &                ln_icedO, rn_hinew, ln_frazil, rn_maxfraz, rn_vfraz, rn_Cfraz, rn_himin,   & 
    545          &                nn_iceflx 
     541      NAMELIST/namice_thd/ ln_icethd, rn_kappa_i, ln_cndi_U64, ln_cndi_P07, ln_dqns_i, rn_cnd_s,   & 
     542         &                 ln_icedH, rn_blow_s,                                                    & 
     543         &                 ln_icedA, rn_beta, rn_dmin,                                             & 
     544         &                 ln_icedO, rn_hinew, ln_frazil, rn_maxfraz, rn_vfraz, rn_Cfraz,          & 
     545         &                 nn_iceflx 
    546546      !!------------------------------------------------------------------- 
    547547      ! 
     
    581581         WRITE(numout,*) '      Threshold relative drift speed for collection of frazil rn_vfraz     = ', rn_vfraz 
    582582         WRITE(numout,*) '      Squeezing coefficient for collection of frazil          rn_Cfraz     = ', rn_Cfraz 
    583          WRITE(numout,*) '   -- iceitd --' 
    584          WRITE(numout,*) '      minimum ice thickness                                   rn_himin     = ', rn_himin  
    585583         WRITE(numout,*) '   -- icestp --' 
    586584         WRITE(numout,*) '      Multicategory heat flux formulation                     nn_iceflx    = ', nn_iceflx 
Note: See TracChangeset for help on using the changeset viewer.