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

Changeset 6746


Ignore:
Timestamp:
2016-06-27T19:20:57+02:00 (8 years ago)
Author:
clem
Message:

landfast ice parameterization + update from trunk + removing useless dom_ice.F90 and limmsh.F90 and limwri_dimg.h90

Location:
branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO
Files:
3 deleted
33 edited

Legend:

Unmodified
Added
Removed

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_2/limmsh_2.F90

    r3625 r6746  
    8080      njeqm1 = njeq - 1  
    8181 
    82       fcor(:,:) = 2. * omega * SIN( gphit(:,:) * rad )   !  coriolis factor at T-point 
     82      fcor(:,:) = 2. * omega * SIN( gphif(:,:) * rad )   !  coriolis factor at T-point 
    8383  
    8484!i    DO jj = 1, jpj 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r6629 r6746  
    203203   INTEGER , PUBLIC ::   nn_icestr        !: ice strength parameterization (0=Hibler79 1=Rothrock75) 
    204204   REAL(wp), PUBLIC ::   rn_pe_rdg        !: ridging work divided by pot. energy change in ridging, nn_icestr = 1 
    205    REAL(wp), PUBLIC ::   rn_pstar         !: determines ice strength (N/M), Hibler JPO79 
     205   REAL(wp), PUBLIC ::   rn_pstar         !: determines ice strength, Hibler JPO79 
    206206   REAL(wp), PUBLIC ::   rn_crhg          !: determines changes in ice strength 
    207207   LOGICAL , PUBLIC ::   ln_icestr_bvf    !: use brine volume to diminish ice strength 
    … …  
    212212   INTEGER , PUBLIC ::   nn_nevp          !: number of iterations for subcycling 
    213213   REAL(wp), PUBLIC ::   rn_relast        !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp)  
     214   LOGICAL , PUBLIC ::   ln_landfast      !: landfast ice parameterization (T or F)  
     215   REAL(wp), PUBLIC ::   rn_gamma         !: fraction of ocean depth that ice must reach to initiate landfast ice 
     216   REAL(wp), PUBLIC ::   rn_icebfr        !: maximum bottom stress per unit area of contact (landfast ice)  
    214217 
    215218   !                                     !!** ice-diffusion namelist (namicehdf) ** 
    … …  
    296299   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   shear_i     !: Shear of the velocity field [s-1] 
    297300   ! 
    298    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sist        !: Average Sea-Ice Surface Temperature [Kelvin] 
    299301   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_bo        !: Sea-Ice bottom temperature [Kelvin]      
    300302   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   frld        !: Leads fraction = 1 - ice fraction 
    … …  
    375377   !                                                                    !  this is an extensive variable that has to be transported 
    376378   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   o_i       !: Sea-Ice Age (days) 
    377    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ov_i      !: Sea-Ice Age times volume per area (days.m) 
    378379   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   oa_i      !: Sea-Ice Age times ice area (days) 
    379380   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   bv_i      !: brine volume 
    … …  
    381382   !! Variables summed over all categories, or associated to all the ice in a single grid cell 
    382383   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   u_ice, v_ice !: components of the ice velocity (m/s) 
    383    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tio_u, tio_v !: components of the ice-ocean stress (N/m2) 
    384384   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vt_i , vt_s  !: ice and snow total volume per unit area (m) 
    385385   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   at_i         !: ice total fractional area (ice concentration) 
    … …  
    393393   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htm_s        !: mean snow thickness over all categories 
    394394   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   om_i         !: mean ice age over all categories 
     395   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tau_icebfr   !: ice friction with bathy (landfast param activated) 
    395396 
    396397   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   t_s      !: Snow temperatures [K] 
    … …  
    466467      ! stay within Fortran's max-line length limit. 
    467468      ii = 1 
    468       ALLOCATE( u_oce    (jpi,jpj) , v_oce    (jpi,jpj) , ahiu      (jpi,jpj) , ahiv     (jpi,jpj) ,  & 
    469          &      hicol    (jpi,jpj) , strength (jpi,jpj) ,                                             & 
    470          &      stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) ,                       & 
    471          &      delta_i  (jpi,jpj) , divu_i   (jpi,jpj) , shear_i   (jpi,jpj) , STAT=ierr(ii) ) 
    472  
    473       ii = ii + 1 
    474       ALLOCATE( sist   (jpi,jpj) , t_bo   (jpi,jpj) ,                                           & 
    475          &      frld   (jpi,jpj) , pfrld  (jpi,jpj) , phicif (jpi,jpj) ,                        & 
     469      ALLOCATE( u_oce   (jpi,jpj) , v_oce    (jpi,jpj) ,                                             & 
     470         &      ahiu    (jpi,jpj) , ahiv     (jpi,jpj) , hicol    (jpi,jpj) ,                        & 
     471         &      strength(jpi,jpj) , stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) ,  & 
     472         &      delta_i (jpi,jpj) , divu_i   (jpi,jpj) , shear_i  (jpi,jpj) , STAT=ierr(ii) ) 
     473 
     474      ii = ii + 1 
     475      ALLOCATE( t_bo   (jpi,jpj) , frld   (jpi,jpj) , pfrld  (jpi,jpj) , phicif (jpi,jpj) ,     & 
    476476         &      wfx_snw(jpi,jpj) , wfx_ice(jpi,jpj) , wfx_sub(jpi,jpj) , wfx_lam(jpi,jpj) ,     & 
    477477         &      wfx_bog(jpi,jpj) , wfx_dyn(jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) ,     & 
    … …  
    493493         &      v_s    (jpi,jpj,jpl) , ht_s  (jpi,jpj,jpl) , t_su  (jpi,jpj,jpl) ,     & 
    494494         &      sm_i   (jpi,jpj,jpl) , smv_i (jpi,jpj,jpl) , o_i   (jpi,jpj,jpl) ,     & 
    495          &      ov_i   (jpi,jpj,jpl) , oa_i  (jpi,jpj,jpl) , bv_i  (jpi,jpj,jpl) ,  STAT=ierr(ii) ) 
    496       ii = ii + 1 
    497       ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) , tio_u(jpi,jpj) , tio_v(jpi,jpj) ,     & 
     495         &      oa_i   (jpi,jpj,jpl) , bv_i  (jpi,jpj,jpl) ,  STAT=ierr(ii) ) 
     496      ii = ii + 1 
     497      ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) ,                                       & 
    498498         &      vt_i (jpi,jpj) , vt_s (jpi,jpj) , at_i (jpi,jpj) , ato_i(jpi,jpj) ,     & 
    499499         &      et_i (jpi,jpj) , et_s (jpi,jpj) , tm_i (jpi,jpj) , bvm_i(jpi,jpj) ,     & 
    500          &      smt_i(jpi,jpj) , tm_su(jpi,jpj) , htm_i(jpi,jpj) , htm_s(jpi,jpj) , om_i(jpi,jpj) , STAT=ierr(ii) ) 
     500         &      smt_i(jpi,jpj) , tm_su(jpi,jpj) , htm_i(jpi,jpj) , htm_s(jpi,jpj) ,     & 
     501         &      om_i (jpi,jpj) , tau_icebfr(jpi,jpj)                              , STAT=ierr(ii) ) 
    501502      ii = ii + 1 
    502503      ALLOCATE( t_s(jpi,jpj,nlay_s,jpl) , e_s(jpi,jpj,nlay_s,jpl) , STAT=ierr(ii) ) 
    … …  
    540541 
    541542      ice_alloc = MAXVAL( ierr(:) ) 
    542       IF( ice_alloc /= 0 )   CALL ctl_warn('ice_alloc_2: failed to allocate arrays.') 
     543      IF( ice_alloc /= 0 )   CALL ctl_warn('ice_alloc: failed to allocate arrays.') 
    543544      ! 
    544545   END FUNCTION ice_alloc 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90

    r5429 r6746  
    1616   !!---------------------------------------------------------------------- 
    1717   USE dom_oce          ! ocean domain 
    18    USE dom_ice          ! LIM-3 domain 
    1918   USE ice              ! LIM-3 variables 
    2019   USE lbclnk           ! lateral boundary condition - MPP exchanges 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limadv_umx.F90

    r6515 r6746  
    1616   USE dom_oce        ! ocean domain 
    1717   USE sbc_oce        ! ocean surface boundary condition 
    18    USE dom_ice        ! ice domain 
    1918   USE ice            ! ice variables 
    2019   ! 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90

    r6515 r6746  
    1818   USE phycst         ! physical constants 
    1919   USE ice            ! LIM-3 variables 
    20    USE dom_ice        ! LIM-3 domain 
    2120   USE dom_oce        ! ocean domain 
    2221   USE in_out_manager ! I/O manager 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90

    r5167 r6746  
    1717   USE ice             ! LIM-3: ice variables 
    1818   USE thd_ice         ! LIM-3: thermodynamical variables 
    19    USE dom_ice         ! LIM-3: ice domain 
    2019   USE sbc_oce         ! Surface boundary condition: ocean fields 
    2120   USE sbc_ice         ! Surface boundary condition: ice   fields 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90

    r6515 r6746  
    1414   !!---------------------------------------------------------------------- 
    1515   USE ice             ! LIM-3: sea-ice variable 
    16    USE dom_ice         ! LIM-3: sea-ice domain 
    1716   USE dom_oce         ! ocean domain 
    1817   USE sbc_oce         ! surface boundary condition: ocean fields 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90

    r6515 r6746  
    1919   USE sbc_ice          ! Surface boundary condition: ice   fields 
    2020   USE ice              ! LIM-3 variables 
    21    USE dom_ice          ! LIM-3 domain 
    2221   USE limrhg           ! LIM-3 rheology 
    2322   USE lbclnk           ! lateral boundary conditions - MPP exchanges 
    … …  
    5857      INTEGER, INTENT(in) ::   kt     ! number of iteration 
    5958      !! 
    60       INTEGER  ::   ji, jj, jl, ja    ! dummy loop indices 
    61       INTEGER  ::   i_j1, i_jpj       ! Starting/ending j-indices for rheology 
    62       REAL(wp), POINTER, DIMENSION(:)   ::   zswitch        ! i-averaged indicator of sea-ice 
    63       REAL(wp), POINTER, DIMENSION(:)   ::   zmsk           ! i-averaged of tmask 
    64       ! 
     59      INTEGER  :: jl, jk ! dummy loop indices 
    6560      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    6661     !!--------------------------------------------------------------------- 
    … …  
    6863      IF( nn_timing == 1 )  CALL timing_start('limdyn') 
    6964 
    70       CALL wrk_alloc( jpj, zswitch, zmsk ) 
    71  
    72       CALL lim_var_agg(1)             ! aggregate ice categories 
     65      CALL lim_var_agg(1)                      ! aggregate ice categories 
    7366 
    7467      IF( kt == nit000 )   CALL lim_dyn_init   ! Initialization (first time-step only) 
    … …  
    7770      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    7871       
     72      ! ice velocities before rheology 
    7973      u_ice_b(:,:) = u_ice(:,:) * umask(:,:,1) 
    8074      v_ice_b(:,:) = v_ice(:,:) * vmask(:,:,1) 
    8175       
     76      ! Landfast ice parameterization: define max bottom friction 
     77      tau_icebfr(:,:) = 0._wp 
     78      IF( ln_landfast ) THEN 
     79         DO jl = 1, jpl 
     80            WHERE( ht_i(:,:,jl) > bathy(:,:) * rn_gamma )  tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 
     81         END DO 
     82      ENDIF 
     83       
    8284      ! Rheology (ice dynamics) 
    83       ! ======== 
    84        
    85       !  Define the j-limits where ice rheology is computed 
    86       ! --------------------------------------------------- 
    87        
    88       IF( lk_mpp .OR. lk_mpp_rep ) THEN                    ! mpp: compute over the whole domain 
    89          i_j1 = 1 
    90          i_jpj = jpj 
    91          IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj ) 
    92          CALL lim_rhg( i_j1, i_jpj ) 
    93       ELSE                                 ! optimization of the computational area 
    94          ! 
    95          DO jj = 1, jpj 
    96             zswitch(jj) = SUM( 1.0 - at_i(:,jj) )   ! = REAL(jpj) if ocean everywhere on a j-line 
    97             zmsk   (jj) = SUM( tmask(:,jj,1)    )   ! = 0         if land  everywhere on a j-line 
    98          END DO 
    99           
    100          IF( l_jeq ) THEN                     ! local domain include both hemisphere 
    101             !                                 ! Rheology is computed in each hemisphere 
    102             !                                 ! only over the ice cover latitude strip 
    103             ! Northern hemisphere 
    104             i_j1  = njeq 
    105             i_jpj = jpj 
    106             DO WHILE ( i_j1 <= jpj .AND. zswitch(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 ) 
    107                i_j1 = i_j1 + 1 
    108             END DO 
    109             i_j1 = MAX( 1, i_j1-2 ) 
    110             IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : NH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj ) 
    111             CALL lim_rhg( i_j1, i_jpj ) 
    112             ! 
    113             ! Southern hemisphere 
    114             i_j1  =  1 
    115             i_jpj = njeq 
    116             DO WHILE ( i_jpj >= 1 .AND. zswitch(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 ) 
    117                i_jpj = i_jpj - 1 
    118             END DO 
    119             i_jpj = MIN( jpj, i_jpj+1 ) 
    120             IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : SH  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj ) 
    121             ! 
    122             CALL lim_rhg( i_j1, i_jpj ) 
    123             ! 
    124          ELSE                                 ! local domain extends over one hemisphere only 
    125             !                                 ! Rheology is computed only over the ice cover 
    126             !                                 ! latitude strip 
    127             i_j1  = 1 
    128             DO WHILE ( i_j1 <= jpj .AND. zswitch(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 ) 
    129                i_j1 = i_j1 + 1 
    130             END DO 
    131             i_j1 = MAX( 1, i_j1-2 ) 
    132              
    133             i_jpj  = jpj 
    134             DO WHILE ( i_jpj >= 1  .AND. zswitch(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 ) 
    135                i_jpj = i_jpj - 1 
    136             END DO 
    137             i_jpj = MIN( jpj, i_jpj+1) 
    138             ! 
    139             IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn  : one hemisphere:  i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj ) 
    140             ! 
    141             CALL lim_rhg( i_j1, i_jpj ) 
    142             ! 
    143          ENDIF 
    144          ! 
    145       ENDIF 
     85      ! ========      
     86      CALL lim_rhg 
    14687      ! 
    147       IF(ln_ctl) THEN   ! Control print 
     88      ! 
     89      ! Control prints 
     90      IF(ln_ctl) THEN 
    14891         CALL prt_ctl_info(' ') 
    14992         CALL prt_ctl_info(' - Cell values : ') 
    … …  
    173116            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_dyn  : sm_i     : ') 
    174117            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_dyn  : smv_i    : ') 
    175             DO ja = 1, nlay_i 
     118            DO jk = 1, nlay_i 
    176119               CALL prt_ctl_info(' ') 
    177                CALL prt_ctl_info(' - Layer : ', ivar1=ja) 
     120               CALL prt_ctl_info(' - Layer : ', ivar1=jk) 
    178121               CALL prt_ctl_info('   ~~~~~~~') 
    179                CALL prt_ctl(tab2d_1=t_i(:,:,ja,jl) , clinfo1= ' lim_dyn  : t_i      : ') 
    180                CALL prt_ctl(tab2d_1=e_i(:,:,ja,jl) , clinfo1= ' lim_dyn  : e_i      : ') 
     122               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_dyn  : t_i      : ') 
     123               CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_dyn  : e_i      : ') 
    181124            END DO 
    182125         END DO 
    … …  
    185128      ! conservation test 
    186129      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    187       ! 
    188       CALL wrk_dealloc( jpj, zswitch, zmsk ) 
    189130      ! 
    190131      IF( nn_timing == 1 )  CALL timing_stop('limdyn') 
    … …  
    207148      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    208149      NAMELIST/namicedyn/ nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio,  rn_creepl, rn_ecc, & 
    209          &                nn_nevp, rn_relast 
     150         &                nn_nevp, rn_relast, ln_landfast, rn_gamma, rn_icebfr 
    210151      !!------------------------------------------------------------------- 
    211152 
    … …  
    223164         WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics ' 
    224165         WRITE(numout,*) '~~~~~~~~~~~~' 
    225          WRITE(numout,*)'    ice strength parameterization (0=Hibler 1=Rothrock)  nn_icestr     = ', nn_icestr  
    226          WRITE(numout,*)'    Including brine volume in ice strength comp.         ln_icestr_bvf = ', ln_icestr_bvf 
    227          WRITE(numout,*)'    Ratio of ridging work to PotEner change in ridging   rn_pe_rdg     = ', rn_pe_rdg  
    228          WRITE(numout,*) '   drag coefficient for oceanic stress                  rn_cio        = ', rn_cio 
    229          WRITE(numout,*) '   first bulk-rheology parameter                        rn_pstar      = ', rn_pstar 
    230          WRITE(numout,*) '   second bulk-rhelogy parameter                        rn_crhg       = ', rn_crhg 
    231          WRITE(numout,*) '   creep limit                                          rn_creepl     = ', rn_creepl 
    232          WRITE(numout,*) '   eccentricity of the elliptical yield curve           rn_ecc        = ', rn_ecc 
    233          WRITE(numout,*) '   number of iterations for subcycling                  nn_nevp       = ', nn_nevp 
    234          WRITE(numout,*) '   ratio of elastic timescale over ice time step        rn_relast     = ', rn_relast 
     166         WRITE(numout,*)'    ice strength parameterization (0=Hibler 1=Rothrock)         nn_icestr     = ', nn_icestr  
     167         WRITE(numout,*)'    Including brine volume in ice strength comp.                ln_icestr_bvf = ', ln_icestr_bvf 
     168         WRITE(numout,*)'    Ratio of ridging work to PotEner change in ridging          rn_pe_rdg     = ', rn_pe_rdg  
     169         WRITE(numout,*) '   drag coefficient for oceanic stress                         rn_cio        = ', rn_cio 
     170         WRITE(numout,*) '   first bulk-rheology parameter                               rn_pstar      = ', rn_pstar 
     171         WRITE(numout,*) '   second bulk-rhelogy parameter                               rn_crhg       = ', rn_crhg 
     172         WRITE(numout,*) '   creep limit                                                 rn_creepl     = ', rn_creepl 
     173         WRITE(numout,*) '   eccentricity of the elliptical yield curve                  rn_ecc        = ', rn_ecc 
     174         WRITE(numout,*) '   number of iterations for subcycling                         nn_nevp       = ', nn_nevp 
     175         WRITE(numout,*) '   ratio of elastic timescale over ice time step               rn_relast     = ', rn_relast 
     176         WRITE(numout,*) '   Landfast: param (T or F)                                    ln_landfast   = ', ln_landfast 
     177         WRITE(numout,*) '   Landfast: fraction of ocean depth that ice must reach       rn_gamma      = ', rn_gamma 
     178         WRITE(numout,*) '   Landfast: maximum bottom stress per unit area of contact    rn_icebfr     = ', rn_icebfr 
    235179      ENDIF 
    236180      ! 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90

    r6693 r6746  
    2323   USE ice              ! sea-ice variables 
    2424   USE par_oce          ! ocean parameters 
    25    USE dom_ice          ! sea-ice domain 
    2625   USE limvar           ! lim_var_salprof 
    2726   USE in_out_manager   ! I/O manager 
    … …  
    153152            DO jj = 1, jpj 
    154153               DO ji = 1, jpi 
    155                   IF( fcor(ji,jj) >= 0._wp ) THEN 
     154                  IF( ff(ji,jj) >= 0._wp ) THEN 
    156155                     zht_i_ini(ji,jj) = rn_hti_ini_n 
    157156                     zht_s_ini(ji,jj) = rn_hts_ini_n 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90

    r6515 r6746  
    1818   USE thd_ice          ! LIM thermodynamics 
    1919   USE ice              ! LIM variables 
    20    USE dom_ice          ! LIM domain 
    2120   USE limvar           ! LIM 
    2221   USE lbclnk           ! lateral boundary condition - MPP exchanges 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90

    r5407 r6746  
    1818   !!   lim_itd_shiftice : 
    1919   !!---------------------------------------------------------------------- 
    20    USE dom_ice          ! LIM-3 domain 
    2120   USE par_oce          ! ocean parameters 
    2221   USE dom_oce          ! ocean domain 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

    r6584 r6746  
    99   !!            3.3  !  2009-05  (G.Garric) addition of the lim2_evp cas 
    1010   !!            3.4  !  2011-01  (A. Porter)  dynamical allocation  
    11    !!            3.5  !  2012-08  (R. Benshila)  AGRIF  
    12    !!---------------------------------------------------------------------- 
    13 #if defined key_lim3 || (  defined key_lim2 && ! defined key_lim2_vp ) 
    14    !!---------------------------------------------------------------------- 
    15    !!   'key_lim3'               OR                     LIM-3 sea-ice model 
    16    !!   'key_lim2' AND NOT 'key_lim2_vp'            EVP LIM-2 sea-ice model 
     11   !!            3.5  !  2012-08  (R. Benshila)  AGRIF 
     12   !!            3.6  !  2016-06  (C. Rousset) Rewriting and add the possibility to use mEVP from Bouillon 2013 
     13   !!---------------------------------------------------------------------- 
     14#if defined key_lim3 
     15   !!---------------------------------------------------------------------- 
     16   !!   'key_lim3'                                      LIM-3 sea-ice model 
    1717   !!---------------------------------------------------------------------- 
    1818   !!   lim_rhg       : computes ice velocities 
    … …  
    2424   USE sbc_oce        ! Surface boundary condition: ocean fields 
    2525   USE sbc_ice        ! Surface boundary condition: ice fields 
    26 #if defined key_lim3 
    27    USE ice            ! LIM-3: ice variables 
    28    USE dom_ice        ! LIM-3: ice domain 
    29    USE limitd_me      ! LIM-3:  
    30 #else 
    31    USE ice_2          ! LIM-2: ice variables 
    32    USE dom_ice_2      ! LIM-2: ice domain 
    33 #endif 
     26   USE ice            ! ice variables 
     27   USE limitd_me      ! ice strength 
    3428   USE lbclnk         ! Lateral Boundary Condition / MPP link 
    3529   USE lib_mpp        ! MPP library 
    … …  
    3832   USE prtctl         ! Print control 
    3933   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    40 #if defined key_agrif && defined key_lim2 
    41    USE agrif_lim2_interp 
    42 #endif 
    43 #if defined key_agrif && defined key_lim3 
     34#if defined key_agrif 
    4435   USE agrif_lim3_interp 
    4536#endif 
    … …  
    5142   PRIVATE 
    5243 
    53    PUBLIC   lim_rhg        ! routine called by lim_dyn (or lim_dyn_2) 
     44   PUBLIC   lim_rhg        ! routine called by lim_dyn 
    5445 
    5546   !! * Substitutions 
    … …  
    6253CONTAINS 
    6354 
    64    SUBROUTINE lim_rhg( k_j1, k_jpj ) 
     55   SUBROUTINE lim_rhg 
    6556      !!------------------------------------------------------------------- 
    6657      !!                 ***  SUBROUTINE lim_rhg  *** 
    … …  
    109100      !!                 e.g. in the Canadian Archipelago 
    110101      !! 
     102      !! ** Notes   : There is the possibility to use mEVP from Bouillon 2013 
     103      !!              (by uncommenting some lines in part 3 and changing alpha and beta parameters) 
     104      !!              but this solution appears very unstable (see Kimmritz et al 2016) 
     105      !! 
    111106      !! References : Hunke and Dukowicz, JPO97 
    112107      !!              Bouillon et al., Ocean Modelling 2009 
     108      !!              Bouillon et al., Ocean Modelling 2013 
    113109      !!------------------------------------------------------------------- 
    114       INTEGER, INTENT(in) ::   k_j1    ! southern j-index for ice computation 
    115       INTEGER, INTENT(in) ::   k_jpj   ! northern j-index for ice computation 
    116       !! 
    117       INTEGER ::   ji, jj   ! dummy loop indices 
    118       INTEGER ::   jter     ! local integers 
     110      INTEGER ::   ji, jj       ! dummy loop indices 
     111      INTEGER ::   jter         ! local integers 
    119112      CHARACTER (len=50) ::   charout 
    120       REAL(wp) ::   zt11, zt12, zt21, zt22, ztagnx, ztagny, delta                         ! 
    121       REAL(wp) ::   za, zstms          ! local scalars 
    122       REAL(wp) ::   zc1, zc2, zc3      ! ice mass 
    123  
    124       REAL(wp) ::   dtevp , z1_dtevp              ! time step for subcycling 
    125       REAL(wp) ::   dtotel, z1_dtotel, ecc2, ecci ! square of yield ellipse eccenticity 
    126       REAL(wp) ::   z0, zr, zcca, zccb            ! temporary scalars 
    127       REAL(wp) ::   zu_ice2, zv_ice1              ! 
    128       REAL(wp) ::   zddc, zdtc                    ! delta on corners and on centre 
    129       REAL(wp) ::   zdst                          ! shear at the center of the grid point 
    130       REAL(wp) ::   zdsshx, zdsshy                ! term for the gradient of ocean surface 
    131       REAL(wp) ::   sigma1, sigma2                ! internal ice stress 
    132  
    133       REAL(wp) ::   zresm         ! Maximal error on ice velocity 
    134       REAL(wp) ::   zintb, zintn  ! dummy argument 
    135  
    136       REAL(wp), POINTER, DIMENSION(:,:) ::   zpresh           ! temporary array for ice strength 
    137       REAL(wp), POINTER, DIMENSION(:,:) ::   zpreshc          ! Ice strength on grid cell corners (zpreshc) 
    138       REAL(wp), POINTER, DIMENSION(:,:) ::   zfrld1, zfrld2   ! lead fraction on U/V points 
    139       REAL(wp), POINTER, DIMENSION(:,:) ::   zmass1, zmass2   ! ice/snow mass on U/V points 
    140       REAL(wp), POINTER, DIMENSION(:,:) ::   zcorl1, zcorl2   ! coriolis parameter on U/V points 
    141       REAL(wp), POINTER, DIMENSION(:,:) ::   za1ct , za2ct    ! temporary arrays 
    142       REAL(wp), POINTER, DIMENSION(:,:) ::   v_oce1           ! ocean u/v component on U points                            
    143       REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce2           ! ocean u/v component on V points 
    144       REAL(wp), POINTER, DIMENSION(:,:) ::   u_ice2, v_ice1   ! ice u/v component on V/U point 
    145       REAL(wp), POINTER, DIMENSION(:,:) ::   zf1   , zf2      ! arrays for internal stresses 
    146       REAL(wp), POINTER, DIMENSION(:,:) ::   zmask            ! mask ocean grid points 
     113 
     114      REAL(wp) ::   dtevp, z1_dtevp                            ! time step for subcycling 
     115      REAL(wp) ::   ecc2, z1_ecc2                              ! square of yield ellipse eccenticity 
     116      REAL(wp) ::   zbeta, zalph1, z1_alph1, zalph2, z1_alph2  ! alpha and beta from Bouillon 2009 and 2013 
     117      REAL(wp) ::   zm1, zm2, zm3                              ! ice/snow mass 
     118      REAL(wp) ::   delta, zp_delf, zds2 
     119      REAL(wp) ::   zTauO, zTauA, zTauB, ZCor, zmt, zu_ice2, zv_ice1, zvel  ! temporary scalars 
     120 
     121      REAL(wp) ::   zsig1, zsig2                               ! internal ice stress 
     122      REAL(wp) ::   zresm                                      ! Maximal error on ice velocity 
     123      REAL(wp) ::   zintb, zintn                               ! dummy argument 
     124 
     125      REAL(wp), POINTER, DIMENSION(:,:) ::   zpresh                    ! temporary array for ice strength 
     126      REAL(wp), POINTER, DIMENSION(:,:) ::   z1_e1t0, z1_e2t0          ! scale factors 
     127      REAL(wp), POINTER, DIMENSION(:,:) ::   zp_delt                   ! P/delta at T points 
     128     ! 
     129      REAL(wp), POINTER, DIMENSION(:,:) ::   za1   , za2               ! ice fraction on U/V points 
     130      REAL(wp), POINTER, DIMENSION(:,:) ::   zmass1, zmass2            ! ice/snow mass on U/V points 
     131      REAL(wp), POINTER, DIMENSION(:,:) ::   zcor1 , zcor2             ! coriolis parameter on U/V points 
     132      REAL(wp), POINTER, DIMENSION(:,:) ::   zspgu , zspgv             ! surface pressure gradient at U/V points 
     133      REAL(wp), POINTER, DIMENSION(:,:) ::   v_oce1, u_oce2, v_ice1, u_ice2  ! ocean/ice u/v component on V/U points                            
     134      REAL(wp), POINTER, DIMENSION(:,:) ::   zf1   , zf2               ! internal stresses 
    147135       
    148       REAL(wp), POINTER, DIMENSION(:,:) ::   zdt              ! tension at centre of grid cells 
    149       REAL(wp), POINTER, DIMENSION(:,:) ::   zds              ! Shear on northeast corner of grid cells 
    150       REAL(wp), POINTER, DIMENSION(:,:) ::   zs1   , zs2      ! Diagonal stress tensor components zs1 and zs2  
    151       REAL(wp), POINTER, DIMENSION(:,:) ::   zs12             ! Non-diagonal stress tensor component zs12 
    152       REAL(wp), POINTER, DIMENSION(:,:) ::   zu_ice, zv_ice, zresr   ! Local error on velocity 
    153       REAL(wp), POINTER, DIMENSION(:,:) ::   zpice            ! array used for the calculation of ice surface slope: 
    154                                                               !   ocean surface (ssh_m) if ice is not embedded 
    155                                                               !   ice top surface if ice is embedded    
    156  
    157       REAL(wp), PARAMETER               ::   zepsi = 1.0e-20_wp ! tolerance parameter 
    158       REAL(wp), PARAMETER               ::   zvmin = 1.0e-03_wp ! ice volume below which ice velocity equals ocean velocity 
     136      REAL(wp), POINTER, DIMENSION(:,:) ::   zdt, zds                  ! tension and shear 
     137      REAL(wp), POINTER, DIMENSION(:,:) ::   zs1, zs2, zs12            ! stress tensor components 
     138      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_ice, zv_ice, zresr     ! check convergence 
     139      REAL(wp), POINTER, DIMENSION(:,:) ::   zpice                     ! array used for the calculation of ice surface slope: 
     140                                                                       !   ocean surface (ssh_m) if ice is not embedded 
     141                                                                       !   ice top surface if ice is embedded    
     142      REAL(wp), POINTER, DIMENSION(:,:) ::   zswitch1, zswitch2        ! dummy arrays 
     143 
     144      REAL(wp), PARAMETER               ::   zepsi = 1.0e-20_wp        ! tolerance parameter 
     145      REAL(wp), PARAMETER               ::   zvmin = 1.0e-03_wp        ! ice volume below which ice velocity equals ocean velocity 
    159146      !!------------------------------------------------------------------- 
    160147 
    161       CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 
    162       CALL wrk_alloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask               ) 
    163       CALL wrk_alloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  ) 
    164       CALL wrk_alloc( jpi,jpj, zs1   , zs2   , zs12   , zresr , zpice                 ) 
    165  
    166 #if  defined key_lim2 && ! defined key_lim2_vp 
    167 # if defined key_agrif 
    168       USE ice_2, vt_s => hsnm 
    169       USE ice_2, vt_i => hicm 
    170 # else 
    171       vt_s => hsnm 
    172       vt_i => hicm 
    173 # endif 
    174       at_i(:,:) = 1. - frld(:,:) 
    175 #endif 
    176 #if defined key_agrif && defined key_lim2  
    177       CALL agrif_rhg_lim2_load      ! First interpolation of coarse values 
    178 #endif 
    179 #if defined key_agrif && defined key_lim3  
     148      CALL wrk_alloc( jpi,jpj, zpresh, z1_e1t0, z1_e2t0, zp_delt ) 
     149      CALL wrk_alloc( jpi,jpj, za1, za2, zmass1, zmass2, zcor1, zcor2 ) 
     150      CALL wrk_alloc( jpi,jpj, zspgu, zspgv, v_oce1, u_oce2, v_ice1, u_ice2, zf1, zf2 ) 
     151      CALL wrk_alloc( jpi,jpj, zds, zdt, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 
     152      CALL wrk_alloc( jpi,jpj, zswitch1, zswitch2 ) 
     153 
     154#if defined key_agrif  
    180155      CALL agrif_interp_lim3('U')   ! First interpolation of coarse values 
    181       CALL agrif_interp_lim3('V')   ! First interpolation of coarse values 
    182 #endif 
    183       ! 
    184       !------------------------------------------------------------------------------! 
    185       ! 1) Ice strength (zpresh)                                ! 
    186       !------------------------------------------------------------------------------! 
    187       ! 
    188       ! Put every vector to 0 
    189       delta_i(:,:) = 0._wp   ; 
    190       zpresh (:,:) = 0._wp   ;   
    191       zpreshc(:,:) = 0._wp 
    192       u_ice2 (:,:) = 0._wp   ;   v_ice1(:,:) = 0._wp 
    193       divu_i (:,:) = 0._wp   ;   zdt   (:,:) = 0._wp   ;   zds(:,:) = 0._wp 
    194       shear_i(:,:) = 0._wp 
    195  
    196 #if defined key_lim3 
    197       CALL lim_itd_me_icestrength( nn_icestr )      ! LIM-3: Ice strength on T-points 
    198 #endif 
    199  
    200       DO jj = k_j1 , k_jpj       ! Ice mass and temp variables 
    201          DO ji = 1 , jpi 
    202 #if defined key_lim3 
    203             zpresh(ji,jj) = tmask(ji,jj,1) *  strength(ji,jj) 
    204 #endif 
    205 #if defined key_lim2 
    206             zpresh(ji,jj) = tmask(ji,jj,1) *  pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) ) 
    207 #endif 
    208             ! zmask = 1 where there is ice or on land 
    209             zmask(ji,jj)    = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - zepsi ) ) ) * tmask(ji,jj,1) 
     156      CALL agrif_interp_lim3('V') 
     157#endif 
     158      ! 
     159      !------------------------------------------------------------------------------! 
     160      ! 0) define some variables 
     161      !------------------------------------------------------------------------------! 
     162      ! ecc2: square of yield ellipse eccenticrity 
     163      ecc2    = rn_ecc * rn_ecc 
     164      z1_ecc2 = 1._wp / ecc2 
     165 
     166      ! Time step for subcycling 
     167      dtevp    = rdt_ice / REAL( nn_nevp ) 
     168      z1_dtevp = 1._wp / dtevp 
     169 
     170      ! alpha parameters (Bouillon 2009) 
     171      zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp 
     172      zalph2 = zalph1 * z1_ecc2 
     173 
     174      ! alpha and beta parameters (Bouillon 2013) 
     175      !!zalph1 = 40. 
     176      !!zalph2 = 40. 
     177      !!zbeta  = 3000. 
     178      !!zbeta = REAL( nn_nevp )   ! close to classical EVP of Hunke (2001) 
     179 
     180      z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 
     181      z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 
     182 
     183      !------------------------------------------------------------------------------! 
     184      ! 1) initialize arrays 
     185      !------------------------------------------------------------------------------! 
     186      ! Initialise stress tensor  
     187      zs1 (:,:) = stress1_i (:,:)  
     188      zs2 (:,:) = stress2_i (:,:) 
     189      zs12(:,:) = stress12_i(:,:) 
     190 
     191      ! Ice strength 
     192      CALL lim_itd_me_icestrength( nn_icestr ) 
     193      zpresh(:,:) = tmask(:,:,1) *  strength(:,:) 
     194 
     195      ! scale factors 
     196      DO jj = 2, jpjm1 
     197         DO ji = fs_2, fs_jpim1 
     198            z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj  ) + e1t(ji,jj  ) ) 
     199            z1_e2t0(ji,jj) = 1._wp / ( e2t(ji  ,jj+1) + e2t(ji,jj  ) ) 
    210200         END DO 
    211201      END DO 
    212  
    213       ! Ice strength on grid cell corners (zpreshc) 
    214       ! needed for calculation of shear stress  
    215       DO jj = k_j1+1, k_jpj-1 
    216          DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1) 
    217             zstms          =  tmask(ji+1,jj+1,1) * wght(ji+1,jj+1,2,2) + tmask(ji,jj+1,1) * wght(ji+1,jj+1,1,2) +   & 
    218                &              tmask(ji+1,jj,1)   * wght(ji+1,jj+1,2,1) + tmask(ji,jj,1)   * wght(ji+1,jj+1,1,1) 
    219             zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) +   & 
    220                &               zpresh(ji+1,jj)   * wght(ji+1,jj+1,2,1) + zpresh(ji,jj)   * wght(ji+1,jj+1,1,1)     & 
    221                &             ) / MAX( zstms, zepsi ) 
    222          END DO 
    223       END DO 
    224       CALL lbc_lnk( zpreshc(:,:), 'F', 1. ) 
     202             
    225203      ! 
    226204      !------------------------------------------------------------------------------! 
    227205      ! 2) Wind / ocean stress, mass terms, coriolis terms 
    228206      !------------------------------------------------------------------------------! 
    229       ! 
    230       !  Wind stress, coriolis and mass terms on the sides of the squares         
    231       !  zfrld1: lead fraction on U-points                                       
    232       !  zfrld2: lead fraction on V-points                                      
    233       !  zmass1: ice/snow mass on U-points                                     
    234       !  zmass2: ice/snow mass on V-points                                    
    235       !  zcorl1: Coriolis parameter on U-points                              
    236       !  zcorl2: Coriolis parameter on V-points                             
    237       !  (ztagnx,ztagny): wind stress on U/V points                        
    238       !  v_oce1: ocean v component on u points                           
    239       !  u_oce2: ocean u component on v points                          
    240207 
    241208      IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==! 
    … …  
    249216         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 
    250217         ! 
    251          zpice(:,:) = ssh_m(:,:) + (  zintn * snwice_mass(:,:) +  zintb * snwice_mass_b(:,:) ) * r1_rau0 
     218         zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0 
    252219         ! 
    253220      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==! 
    … …  
    255222      ENDIF 
    256223 
    257       DO jj = k_j1+1, k_jpj-1 
     224      DO jj = 2, jpjm1 
    258225         DO ji = fs_2, fs_jpim1 
    259226 
    260             zc1 = tmask(ji  ,jj  ,1) * ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) ) 
    261             zc2 = tmask(ji+1,jj  ,1) * ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) ) 
    262             zc3 = tmask(ji  ,jj+1,1) * ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) ) 
    263  
    264             zt11 = tmask(ji  ,jj,1) * e1t(ji  ,jj) 
    265             zt12 = tmask(ji+1,jj,1) * e1t(ji+1,jj) 
    266             zt21 = tmask(ji,jj  ,1) * e2t(ji,jj  ) 
    267             zt22 = tmask(ji,jj+1,1) * e2t(ji,jj+1) 
    268  
    269             ! Leads area. 
    270             zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + zepsi ) 
    271             zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + zepsi ) 
    272  
    273             ! Mass, coriolis coeff. and currents 
    274             zmass1(ji,jj) = ( zt12 * zc1 + zt11 * zc2 ) / ( zt11 + zt12 + zepsi ) 
    275             zmass2(ji,jj) = ( zt22 * zc1 + zt21 * zc3 ) / ( zt21 + zt22 + zepsi ) 
    276             zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj) * fcor(ji,jj) + e1t(ji,jj) * fcor(ji+1,jj) )   & 
    277                &                          / ( e1t(ji,jj) + e1t(ji+1,jj) + zepsi ) 
    278             zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1) * fcor(ji,jj) + e2t(ji,jj) * fcor(ji,jj+1) )   & 
    279                &                          / ( e2t(ji,jj+1) + e2t(ji,jj) + zepsi ) 
    280             ! 
    281             ! Ocean has no slip boundary condition 
    282             v_oce1(ji,jj)  = 0.5 * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji,jj)      & 
    283                &                   + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji+1,jj) )  & 
    284                &                   / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)   
    285  
    286             u_oce2(ji,jj)  = 0.5 * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj)      & 
    287                &                   + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj+1) )  & 
    288                &                   / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 
    289  
    290             ! Wind stress at U,V-point 
    291             ztagnx = ( 1. - zfrld1(ji,jj) ) * utau_ice(ji,jj) 
    292             ztagny = ( 1. - zfrld2(ji,jj) ) * vtau_ice(ji,jj) 
    293  
    294             ! Computation of the velocity field taking into account the ice internal interaction. 
    295             ! Terms that are independent of the velocity field. 
    296  
    297             ! SB On utilise maintenant le gradient de la pente de l'ocean 
    298             ! include it later 
    299  
    300             zdsshx =  ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj) 
    301             zdsshy =  ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj) 
    302  
    303             za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx 
    304             za2ct(ji,jj) = ztagny - zmass2(ji,jj) * grav * zdsshy 
     227            ! ice fraction at U-V points 
     228            za1(ji,jj) = ( at_i(ji,jj) * e1t(ji+1,jj) + at_i(ji+1,jj) * e1t(ji,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
     229            za2(ji,jj) = ( at_i(ji,jj) * e2t(ji,jj+1) + at_i(ji,jj+1) * e2t(ji,jj) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
     230 
     231            ! Ice/snow mass at U-V points 
     232            zm1 = ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) ) 
     233            zm2 = ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) ) 
     234            zm3 = ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) ) 
     235            zmass1(ji,jj) = ( zm1 * e1t(ji+1,jj) + zm2 * e1t(ji,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
     236            zmass2(ji,jj) = ( zm1 * e2t(ji,jj+1) + zm3 * e2t(ji,jj) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
     237 
     238            ! Ocean currents at U-V points 
     239            v_oce1(ji,jj)  = 0.5_wp * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji+1,jj)    & 
     240               &                      + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
     241             
     242            u_oce2(ji,jj)  = 0.5_wp * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj+1)    & 
     243               &                      + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
     244 
     245            ! Coriolis at U-V points (m*f) 
     246            zcor1(ji,jj) =   zmass1(ji,jj) * 0.5_wp * ( ff(ji,jj) + ff(ji  ,jj-1) ) 
     247            zcor2(ji,jj) = - zmass2(ji,jj) * 0.5_wp * ( ff(ji,jj) + ff(ji-1,jj  ) ) 
     248 
     249            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points 
     250            zspgu(ji,jj) = - zmass1(ji,jj) * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj) 
     251            zspgv(ji,jj) = - zmass2(ji,jj) * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj) 
     252 
     253            ! switches 
     254            zswitch1(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) 
     255            zswitch2(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) 
    305256 
    306257         END DO 
    … …  
    312263      !------------------------------------------------------------------------------! 
    313264      ! 
    314       ! Time step for subcycling 
    315       dtevp  = rdt_ice / nn_nevp 
    316 #if defined key_lim3 
    317       dtotel = dtevp / ( 2._wp * rn_relast * rdt_ice ) 
    318 #else 
    319       dtotel = dtevp / ( 2._wp * telast ) 
    320 #endif 
    321       z1_dtotel = 1._wp / ( 1._wp + dtotel ) 
    322       z1_dtevp  = 1._wp / dtevp 
    323       !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter) 
    324       ecc2 = rn_ecc * rn_ecc 
    325       ecci = 1. / ecc2 
    326  
    327       !-Initialise stress tensor  
    328       zs1 (:,:) = stress1_i (:,:)  
    329       zs2 (:,:) = stress2_i (:,:) 
    330       zs12(:,:) = stress12_i(:,:) 
    331  
    332265      !                                               !----------------------! 
    333266      DO jter = 1 , nn_nevp                           !    loop over jter    ! 
    334267         !                                            !----------------------!         
    335          DO jj = k_j1, k_jpj-1 
    336             zu_ice(:,jj) = u_ice(:,jj)    ! velocity at previous time step 
    337             zv_ice(:,jj) = v_ice(:,jj) 
    338          END DO 
    339  
    340          DO jj = k_j1+1, k_jpj-1 
    341             DO ji = fs_2, fs_jpim1   !RB bug no vect opt due to zmask 
    342  
    343                !   
    344                !- Divergence, tension and shear (Section a. Appendix B of Hunke & Dukowicz, 2002) 
    345                !- divu_i(:,:), zdt(:,:): divergence and tension at centre of grid cells 
    346                !- zds(:,:): shear on northeast corner of grid cells 
    347                ! 
    348                !- IMPORTANT REMINDER: Dear Gurvan, note that, the way these terms are coded,  
    349                !                      there are many repeated calculations.  
    350                !                      Speed could be improved by regrouping terms. For 
    351                !                      the moment, however, the stress is on clarity of coding to avoid 
    352                !                      bugs (Martin, for Miguel). 
    353                ! 
    354                !- ALSO: arrays zdt, zds and delta could  
    355                !  be removed in the future to minimise memory demand. 
    356                ! 
    357                !- MORE NOTES: Note that we are calculating deformation rates and stresses on the corners of 
    358                !              grid cells, exactly as in the B grid case. For simplicity, the indexation on 
    359                !              the corners is the same as in the B grid. 
    360                ! 
    361                ! 
    362                divu_i(ji,jj) = (  e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
    363                   &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
     268         ! Convergence test 
     269         IF(ln_ctl) THEN 
     270            DO jj = 1, jpjm1 
     271               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 
     272               zv_ice(:,jj) = v_ice(:,jj) 
     273            END DO 
     274         ENDIF 
     275 
     276         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 
     277         DO jj = 2, jpjm1 
     278            DO ji = fs_2, fs_jpim1 
     279                
     280               ! divergence at T points 
     281               divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     282                  &            + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
    364283                  &            ) * r1_e12t(ji,jj) 
    365284 
     285               ! tension at T points 
    366286               zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
    367287                  &         - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
    368288                  &         ) * r1_e12t(ji,jj) 
    369289 
    370                ! 
     290               ! shear at F points 
    371291               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)   & 
    372292                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
    373                   &         ) * r1_e12f(ji,jj) * ( 2._wp - fmask(ji,jj,1) )   & 
    374                   &         * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1) 
    375  
    376  
    377                v_ice1(ji,jj)  = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     & 
    378                   &                      + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) )   & 
    379                   &                      / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)  
    380  
    381                u_ice2(ji,jj)  = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
    382                   &                      + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) )   & 
    383                   &                      / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 
    384             END DO 
    385          END DO 
    386  
    387          CALL lbc_lnk_multi( v_ice1, 'U', -1., u_ice2, 'V', -1. )      ! lateral boundary cond. 
     293                  &         ) * r1_e12f(ji,jj) 
     294 
     295               ! u_ice at V point 
     296               u_ice2(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
     297                  &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
     298                
     299               ! v_ice at U point 
     300               v_ice1(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     & 
     301                  &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
     302 
     303            END DO 
     304         END DO 
     305         CALL lbc_lnk( zds, 'F', 1. ) 
     306 
     307         DO jj = 2, jpjm1 
     308            DO ji = 2, jpim1 ! no vector loop 
     309 
     310               ! shear**2 at T points (doc eq. A16) 
     311               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e12f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e12f(ji-1,jj  )  & 
     312                  &   + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1)  & 
     313                  &   ) * 0.25_wp * r1_e12t(ji,jj) 
     314               
     315               ! delta at T points 
     316               delta          = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt(ji,jj) * zdt(ji,jj) + zds2 ) * usecc2 )   
     317               delta_i(ji,jj) = delta + rn_creepl 
     318 
     319               ! P/delta at T points 
     320               zp_delt(ji,jj) = zpresh(ji,jj) / delta_i(ji,jj) 
     321            END DO 
     322         END DO 
     323         CALL lbc_lnk( zp_delt, 'T', 1. ) 
     324 
     325         ! --- Stress tensor --- ! 
     326         DO jj = 2, jpjm1 
     327            DO ji = 2, jpim1 ! no vector loop 
     328                
     329               ! P/delta at F points 
     330               zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) ) 
     331                
     332               ! stress tensor at T points 
     333               zs1(ji,jj) = ( zs1 (ji,jj) * zalph1 + zp_delt(ji,jj) * ( divu_i(ji,jj) - (delta_i(ji,jj)-rn_creepl) ) ) * z1_alph1 
     334               zs2(ji,jj) = ( zs2 (ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt(ji,jj) * z1_ecc2  )                      ) * z1_alph2 
     335               ! F points 
     336               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf        * ( zds(ji,jj) * z1_ecc2  ) * 0.5_wp             ) * z1_alph2 
     337            END DO 
     338         END DO 
     339         CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 
     340  
     341         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- ! 
     342         DO jj = 2, jpjm1 
     343            DO ji = fs_2, fs_jpim1                
     344               ! U points 
     345               zf1(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             & 
     346                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    & 
     347                  &                    ) * r1_e2u(ji,jj)                                                                      & 
     348                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  & 
     349                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              & 
     350                  &                  ) * r1_e12u(ji,jj) 
     351 
     352               ! V points 
     353               zf2(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             & 
     354                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    & 
     355                  &                    ) * r1_e1v(ji,jj)                                                                      & 
     356                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  & 
     357                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              & 
     358                  &                  ) * r1_e12v(ji,jj) 
     359            END DO 
     360         END DO 
     361         ! 
     362         ! --- Computation of ice velocity --- ! 
     363         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen smart and large nn_nevp 
     364         !  Bouillon et al. 2009 (eq 34-35) => stable 
     365         IF( MOD(jter,2) .EQ. 0 ) THEN ! even iterations 
     366             
     367            DO jj = 2, jpjm1 
     368               DO ji = fs_2, fs_jpim1 
     369 
     370                  zvel    = SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
     371                     &          + ( u_ice2(ji,jj) - u_oce2(ji,jj) ) * ( u_ice2(ji,jj) - u_oce2(ji,jj) ) ) 
     372 
     373                  zTauA   =   za2(ji,jj) * vtau_ice(ji,jj) 
     374                  zTauO   =   za2(ji,jj) * rhoco * zvel 
     375                  ztauB   = - tau_icebfr(ji,jj) / zvel 
     376                  zCor    =   zcor2(ji,jj)  * u_ice2(ji,jj) 
     377                  zmt     =   zmass2(ji,jj) * z1_dtevp 
     378                   
     379                  ! Bouillon 2009 
     380                  v_ice(ji,jj) = ( zmt * v_ice(ji,jj) + zf2(ji,jj) + zCor + zTauA + zTauO * v_oce(ji,jj) + zspgv(ji,jj)  & 
     381                     &           ) / MAX( zmt + zTauO - zTauB, zepsi ) * zswitch2(ji,jj) 
     382                  ! Bouillon 2013 
     383                  !!v_ice(ji,jj) = ( zmt * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  & 
     384                  !!   &           + zf2(ji,jj) + zCor + zTauA + zTauO * v_oce(ji,jj) + zspgv(ji,jj)  & 
     385                  !!   &           ) / MAX( zmt * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch2(ji,jj) 
     386 
     387               END DO 
     388            END DO 
     389            CALL lbc_lnk( v_ice, 'V', -1. ) 
     390             
     391#if defined key_agrif 
     392            CALL agrif_interp_lim3( 'V' ) 
     393#endif 
     394#if defined key_bdy 
     395            CALL bdy_ice_lim_dyn( 'V' ) 
     396#endif          
     397 
     398            DO jj = 2, jpjm1 
     399               DO ji = fs_2, fs_jpim1 
     400 
     401                  ! v_ice at U point 
     402                  zv_ice1 = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     & 
     403                     &               + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)  
     404                   
     405                  zvel    = SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
     406                     &          + ( v_ice1(ji,jj) - v_oce1(ji,jj) ) * ( v_ice1(ji,jj) - v_oce1(ji,jj) ) ) 
     407 
     408                  zTauA   =   za1(ji,jj) * utau_ice(ji,jj) 
     409                  zTauO   =   za1(ji,jj) * rhoco * zvel 
     410                  ztauB   = - tau_icebfr(ji,jj) / zvel 
     411                  zCor    =   zcor1(ji,jj)  * zv_ice1 
     412                  zmt     =   zmass1(ji,jj) * z1_dtevp 
     413                   
     414                  ! Bouillon 2009 
     415                  u_ice(ji,jj) = ( zmt * u_ice(ji,jj) + zf1(ji,jj) + zCor + zTauA + zTauO * u_oce(ji,jj) + zspgu(ji,jj)  & 
     416                     &           ) / MAX( zmt + zTauO - zTauB, zepsi ) * zswitch1(ji,jj) 
     417                  ! Bouillon 2013 
     418                  !!u_ice(ji,jj) = ( zmt * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  & 
     419                  !!   &           + zf1(ji,jj) + zCor + zTauA + zTauO * u_oce(ji,jj) + zspgu(ji,jj)  & 
     420                  !!   &           ) / MAX( zmt * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch1(ji,jj) 
     421               END DO 
     422            END DO 
     423            CALL lbc_lnk( u_ice, 'U', -1. ) 
     424             
     425#if defined key_agrif 
     426            CALL agrif_interp_lim3( 'U' ) 
     427#endif 
     428#if defined key_bdy 
     429            CALL bdy_ice_lim_dyn( 'U' ) 
     430#endif          
     431 
     432         ELSE ! odd iterations 
     433 
     434            DO jj = 2, jpjm1 
     435               DO ji = fs_2, fs_jpim1 
     436 
     437                  zvel    = SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
     438                     &          + ( v_ice1(ji,jj) - v_oce1(ji,jj) ) * ( v_ice1(ji,jj) - v_oce1(ji,jj) ) ) 
     439 
     440                  zTauA   =   za1(ji,jj) * utau_ice(ji,jj) 
     441                  zTauO   =   za1(ji,jj) * rhoco * zvel 
     442                  ztauB   = - tau_icebfr(ji,jj) / zvel 
     443                  zCor    =   zcor1(ji,jj)  * v_ice1(ji,jj) 
     444                  zmt     =   zmass1(ji,jj) * z1_dtevp 
     445 
     446                  ! Bouillon 2009 
     447                  u_ice(ji,jj) = ( zmt * u_ice(ji,jj) + zf1(ji,jj) + zCor + zTauA + zTauO * u_oce(ji,jj) + zspgu(ji,jj)  & 
     448                     &           ) / MAX( zmt + zTauO - zTauB, zepsi ) * zswitch1(ji,jj) 
     449                  ! Bouillon 2013 
     450                  !!u_ice(ji,jj) = ( zmt * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  & 
     451                  !!   &           + zf1(ji,jj) + zCor + zTauA + zTauO * u_oce(ji,jj) + zspgu(ji,jj)  & 
     452                  !!   &           ) / MAX( zmt * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch1(ji,jj) 
     453               END DO 
     454            END DO 
     455            CALL lbc_lnk( u_ice, 'U', -1. ) 
     456             
     457#if defined key_agrif 
     458            CALL agrif_interp_lim3( 'U' ) 
     459#endif 
     460#if defined key_bdy 
     461            CALL bdy_ice_lim_dyn( 'U' ) 
     462#endif          
     463 
     464            DO jj = 2, jpjm1 
     465               DO ji = fs_2, fs_jpim1 
     466 
     467                  ! u_ice at V point 
     468                  zu_ice2 = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
     469                     &               + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
     470 
     471                  zvel    = SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
     472                     &          + ( u_ice2(ji,jj) - u_oce2(ji,jj) ) * ( u_ice2(ji,jj) - u_oce2(ji,jj) ) ) 
    388473          
    389          DO jj = k_j1+1, k_jpj-1 
    390             DO ji = fs_2, fs_jpim1 
    391  
    392                !- Calculate Delta at centre of grid cells 
    393                zdst          = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )   & 
    394                   &            + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1)   & 
    395                   &            ) * r1_e12t(ji,jj) 
    396  
    397                delta          = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 )   
    398                delta_i(ji,jj) = delta + rn_creepl 
    399  
    400                !- Calculate Delta on corners 
    401                zddc  = (  ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  & 
    402                   &     + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  & 
    403                   &    ) * r1_e12f(ji,jj) 
    404  
    405                zdtc  = (- ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  & 
    406                   &     + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  & 
    407                   &    ) * r1_e12f(ji,jj) 
    408  
    409                zddc = SQRT( zddc**2 + ( zdtc**2 + zds(ji,jj)**2 ) * usecc2 ) + rn_creepl 
    410  
    411                !-Calculate stress tensor components zs1 and zs2 at centre of grid cells (see section 3.5 of CICE user's guide). 
    412                zs1(ji,jj)  = ( zs1 (ji,jj) + dtotel * ( divu_i(ji,jj) - delta ) / delta_i(ji,jj)   * zpresh(ji,jj)   & 
    413                   &          ) * z1_dtotel 
    414                zs2(ji,jj)  = ( zs2 (ji,jj) + dtotel *         ecci * zdt(ji,jj) / delta_i(ji,jj)   * zpresh(ji,jj)   & 
    415                   &          ) * z1_dtotel 
    416                !-Calculate stress tensor component zs12 at corners 
    417                zs12(ji,jj) = ( zs12(ji,jj) + dtotel *         ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj)  & 
    418                   &          ) * z1_dtotel  
    419  
    420             END DO 
    421          END DO 
    422  
    423          CALL lbc_lnk_multi( zs1 , 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 
    424   
    425          ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) 
    426          DO jj = k_j1+1, k_jpj-1 
    427             DO ji = fs_2, fs_jpim1 
    428                !- contribution of zs1, zs2 and zs12 to zf1 
    429                zf1(ji,jj) = 0.5 * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)  & 
    430                   &             + ( zs2(ji+1,jj) * e2t(ji+1,jj)**2 - zs2(ji,jj) * e2t(ji,jj)**2 ) * r1_e2u(ji,jj)          & 
    431                   &             + 2.0 * ( zs12(ji,jj) * e1f(ji,jj)**2 - zs12(ji,jj-1) * e1f(ji,jj-1)**2 ) * r1_e1u(ji,jj)  & 
    432                   &                ) * r1_e12u(ji,jj) 
    433                ! contribution of zs1, zs2 and zs12 to zf2 
    434                zf2(ji,jj) = 0.5 * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)  & 
    435                   &             - ( zs2(ji,jj+1) * e1t(ji,jj+1)**2 - zs2(ji,jj) * e1t(ji,jj)**2 ) * r1_e1v(ji,jj)          & 
    436                   &             + 2.0 * ( zs12(ji,jj) * e2f(ji,jj)**2 - zs12(ji-1,jj) * e2f(ji-1,jj)**2 ) * r1_e2v(ji,jj)  & 
    437                   &               )  * r1_e12v(ji,jj) 
    438             END DO 
    439          END DO 
    440          ! 
    441          ! Computation of ice velocity 
    442          ! 
    443          ! Both the Coriolis term and the ice-ocean drag are solved semi-implicitly. 
    444          ! 
    445          IF (MOD(jter,2).eq.0) THEN  
    446  
    447             DO jj = k_j1+1, k_jpj-1 
    448                DO ji = fs_2, fs_jpim1 
    449                   rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1) 
    450                   z0           = zmass1(ji,jj) * z1_dtevp 
    451  
    452                   ! SB modif because ocean has no slip boundary condition 
    453                   zv_ice1      = 0.5 * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji  ,jj)     & 
    454                      &                 + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) )   & 
    455                      &                 / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 
    456                   za           = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 +  & 
    457                      &                         ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) ) 
    458                   zr           = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj) 
    459                   zcca         = z0 + za 
    460                   zccb         = zcorl1(ji,jj) 
    461                   u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch  
     474                  zTauA   =   za2(ji,jj) * vtau_ice(ji,jj) 
     475                  zTauO   =   za2(ji,jj) * rhoco * zvel 
     476                  ztauB   = - tau_icebfr(ji,jj) / zvel 
     477                  zCor    =   zcor2(ji,jj)  * zu_ice2 
     478                  zmt     =   zmass2(ji,jj) * z1_dtevp 
     479                   
     480                  ! Bouillon 2009 
     481                  v_ice(ji,jj) = ( zmt * v_ice(ji,jj) + zf2(ji,jj) + zCor + zTauA + zTauO * v_oce(ji,jj) + zspgv(ji,jj)  & 
     482                     &           ) / MAX( zmt + zTauO - zTauB, zepsi ) * zswitch2(ji,jj) 
     483                  ! Bouillon 2013 
     484                  !!v_ice(ji,jj) = ( zmt * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  & 
     485                  !!   &           + zf2(ji,jj) + zCor + zTauA + zTauO * v_oce(ji,jj) + zspgv(ji,jj)  & 
     486                  !!   &           ) / MAX( zmt * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch2(ji,jj) 
     487 
    462488               END DO 
    463489            END DO 
    464  
    465             CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 
    466 #if defined key_agrif && defined key_lim2 
    467             CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 
    468 #endif 
    469 #if defined key_agrif && defined key_lim3 
    470             CALL agrif_interp_lim3( 'U' ) 
     490            CALL lbc_lnk( v_ice, 'V', -1. ) 
     491             
     492#if defined key_agrif 
     493            CALL agrif_interp_lim3( 'V' ) 
    471494#endif 
    472495#if defined key_bdy 
    473          CALL bdy_ice_lim_dyn( 'U' ) 
    474 #endif          
    475  
    476             DO jj = k_j1+1, k_jpj-1 
    477                DO ji = fs_2, fs_jpim1 
    478  
    479                   rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1) 
    480                   z0           = zmass2(ji,jj) * z1_dtevp 
    481                   ! SB modif because ocean has no slip boundary condition 
    482                   zu_ice2      = 0.5 * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj)       & 
    483                      &                 + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) )   & 
    484                      &                 / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 
    485                   za           = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 +  &  
    486                      &                         ( v_ice(ji,jj) - v_oce(ji,jj))**2 ) * ( 1.0 - zfrld2(ji,jj) ) 
    487                   zr           = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj) 
    488                   zcca         = z0 + za 
    489                   zccb         = zcorl2(ji,jj) 
    490                   v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch 
    491                END DO 
    492             END DO 
    493  
    494             CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 
    495 #if defined key_agrif && defined key_lim2 
    496             CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 
    497 #endif 
    498 #if defined key_agrif && defined key_lim3 
    499             CALL agrif_interp_lim3( 'V' ) 
    500 #endif 
    501 #if defined key_bdy 
    502          CALL bdy_ice_lim_dyn( 'V' ) 
    503 #endif          
    504  
    505          ELSE  
    506             DO jj = k_j1+1, k_jpj-1 
    507                DO ji = fs_2, fs_jpim1 
    508                   rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1) 
    509                   z0           = zmass2(ji,jj) * z1_dtevp 
    510                   ! SB modif because ocean has no slip boundary condition 
    511                   zu_ice2      = 0.5 * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj)       & 
    512                      &                  +( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) )   & 
    513                      &                 / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)    
    514  
    515                   za           = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 +  & 
    516                      &                         ( v_ice(ji,jj) - v_oce(ji,jj) )**2 ) * ( 1.0 - zfrld2(ji,jj) ) 
    517                   zr           = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj) 
    518                   zcca         = z0 + za 
    519                   zccb         = zcorl2(ji,jj) 
    520                   v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch 
    521                END DO 
    522             END DO 
    523  
    524             CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 
    525 #if defined key_agrif && defined key_lim2 
    526             CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 
    527 #endif 
    528 #if defined key_agrif && defined key_lim3 
    529             CALL agrif_interp_lim3( 'V' ) 
    530 #endif 
    531 #if defined key_bdy 
    532          CALL bdy_ice_lim_dyn( 'V' ) 
    533 #endif          
    534  
    535             DO jj = k_j1+1, k_jpj-1 
    536                DO ji = fs_2, fs_jpim1 
    537                   rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1) 
    538                   z0           = zmass1(ji,jj) * z1_dtevp 
    539                   zv_ice1      = 0.5 * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji,jj)       & 
    540                      &                 + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) )   & 
    541                      &                 / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 
    542  
    543                   za           = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 +  & 
    544                      &                         ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) ) 
    545                   zr           = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj) 
    546                   zcca         = z0 + za 
    547                   zccb         = zcorl1(ji,jj) 
    548                   u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch  
    549                END DO 
    550             END DO 
    551  
    552             CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 
    553 #if defined key_agrif && defined key_lim2 
    554             CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 
    555 #endif 
    556 #if defined key_agrif && defined key_lim3 
    557             CALL agrif_interp_lim3( 'U' ) 
    558 #endif 
    559 #if defined key_bdy 
    560          CALL bdy_ice_lim_dyn( 'U' ) 
     496            CALL bdy_ice_lim_dyn( 'V' ) 
    561497#endif          
    562498 
    563499         ENDIF 
    564500          
     501         ! Convergence test 
    565502         IF(ln_ctl) THEN 
    566             !---  Convergence test. 
    567             DO jj = k_j1+1 , k_jpj-1 
     503            DO jj = 2 , jpjm1 
    568504               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 
    569505            END DO 
    570             zresm = MAXVAL( zresr( 1:jpi, k_j1+1:k_jpj-1 ) ) 
     506            zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) ) 
    571507            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain 
    572508         ENDIF 
    573  
     509         ! 
    574510         !                                                ! ==================== ! 
    575511      END DO                                              !  end loop over jter  ! 
    … …  
    577513      ! 
    578514      !------------------------------------------------------------------------------! 
    579       ! 4) Prevent ice velocities when the ice is thin 
     515      ! 4) Limit ice velocities when ice is thin 
    580516      !------------------------------------------------------------------------------! 
    581517      ! If the ice volume is below zvmin then ice velocity should equal the 
    582518      ! ocean velocity. This prevents high velocity when ice is thin 
    583       DO jj = k_j1+1, k_jpj-1 
     519       DO jj = 2, jpjm1 
    584520         DO ji = fs_2, fs_jpim1 
    585             IF ( vt_i(ji,jj) <= zvmin ) THEN 
    586                u_ice(ji,jj) = u_oce(ji,jj) 
    587                v_ice(ji,jj) = v_oce(ji,jj) 
    588             ENDIF 
     521            rswitch      = MAX( 0._wp, SIGN( 1._wp, vt_i(ji,jj) - zvmin ) ) 
     522            u_ice(ji,jj) = ( rswitch * u_ice(ji,jj) + (1._wp - rswitch) * u_oce(ji,jj) ) * zswitch1(ji,jj) 
     523            v_ice(ji,jj) = ( rswitch * v_ice(ji,jj) + (1._wp - rswitch) * v_oce(ji,jj) ) * zswitch2(ji,jj) 
    589524         END DO 
    590525      END DO 
    591526 
    592       CALL lbc_lnk_multi( u_ice(:,:), 'U', -1., v_ice(:,:), 'V', -1. ) 
    593  
    594 #if defined key_agrif && defined key_lim2 
    595       CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'U' ) 
    596       CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'V' ) 
    597 #endif 
    598 #if defined key_agrif && defined key_lim3 
     527      CALL lbc_lnk_multi( u_ice, 'U', -1., v_ice, 'V', -1. ) 
     528 
     529#if defined key_agrif 
    599530      CALL agrif_interp_lim3( 'U' ) 
    600531      CALL agrif_interp_lim3( 'V' ) 
    … …  
    605536#endif          
    606537 
    607       DO jj = k_j1+1, k_jpj-1  
     538      !------------------------------------------------------------------------------! 
     539      ! 5) Recompute delta, shear and div (inputs for mechanical redistribution)  
     540      !------------------------------------------------------------------------------! 
     541 
     542      ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 
     543      DO jj = 2, jpjm1 
    608544         DO ji = fs_2, fs_jpim1 
    609             IF ( vt_i(ji,jj) <= zvmin ) THEN 
    610                v_ice1(ji,jj)  = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji,  jj-1) ) * e1t(ji+1,jj)     & 
    611                   &                      + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) )   & 
    612                   &                      / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 
    613  
    614                u_ice2(ji,jj)  = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
    615                   &                      + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) )   & 
    616                   &                      / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 
    617             ENDIF  
     545             
     546            ! divergence at T points 
     547            divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     548               &            + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
     549               &            ) * r1_e12t(ji,jj) 
     550             
     551            ! tension at T points 
     552            zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
     553               &         - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
     554               &         ) * r1_e12t(ji,jj) 
     555             
     556            ! shear at F points 
     557            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)   & 
     558               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
     559               &         ) * r1_e12f(ji,jj) 
     560             
    618561         END DO 
    619562      END DO 
    620  
    621       CALL lbc_lnk_multi( u_ice2(:,:), 'V', -1., v_ice1(:,:), 'U', -1. ) 
    622  
    623       ! Recompute delta, shear and div, inputs for mechanical redistribution  
    624       DO jj = k_j1+1, k_jpj-1 
    625          DO ji = fs_2, jpim1   !RB bug no vect opt due to zmask 
    626             !- divu_i(:,:), zdt(:,:): divergence and tension at centre  
    627             !- zds(:,:): shear on northeast corner of grid cells 
    628             IF ( vt_i(ji,jj) <= zvmin ) THEN 
    629  
    630                divu_i(ji,jj) = (  e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj  ) * u_ice(ji-1,jj  )   & 
    631                   &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji  ,jj-1) * v_ice(ji  ,jj-1)   & 
    632                   &            ) * r1_e12t(ji,jj) 
    633  
    634                zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)  & 
    635                   &          -( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)  & 
    636                   &         ) * r1_e12t(ji,jj) 
    637                ! 
    638                ! SB modif because ocean has no slip boundary condition  
    639                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)  & 
    640                   &          +( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  & 
    641                   &         ) * r1_e12f(ji,jj) * ( 2.0 - fmask(ji,jj,1) )                                     & 
    642                   &         * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1) 
    643  
    644                zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )    & 
    645                   &   + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1) ) * r1_e12t(ji,jj) 
    646  
    647                delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 )   
    648                delta_i(ji,jj) = delta + rn_creepl 
    649              
    650             ENDIF 
     563      CALL lbc_lnk_multi( divu_i, 'T', 1., zds, 'F', 1. ) 
     564       
     565 
     566      DO jj = 2, jpjm1 
     567         DO ji = 2, jpim1 ! no vector loop 
     568             
     569            ! shear**2 at T points (doc eq. A16) 
     570            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e12f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e12f(ji-1,jj  )  & 
     571               &   + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1)  & 
     572               &   ) * 0.25_wp * r1_e12t(ji,jj) 
     573             
     574            ! delta at T points 
     575            delta          = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zds2 ) * usecc2 )   
     576            delta_i(ji,jj) = delta + rn_creepl 
     577             
     578            shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zds2 ) 
    651579         END DO 
    652580      END DO 
    653       ! 
    654       !------------------------------------------------------------------------------! 
    655       ! 5) Store stress tensor and its invariants 
    656       !------------------------------------------------------------------------------! 
    657       ! * Invariants of the stress tensor are required for limitd_me 
    658       !   (accelerates convergence and improves stability) 
    659       DO jj = k_j1+1, k_jpj-1 
    660          DO ji = fs_2, fs_jpim1 
    661             zdst           = (  e2u(ji,jj) * v_ice1(ji,jj) - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)  &    
    662                &              + e1v(ji,jj) * u_ice2(ji,jj) - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1) ) * r1_e12t(ji,jj)  
    663             shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst ) 
    664          END DO 
    665       END DO 
    666  
    667       ! Lateral boundary condition 
    668       CALL lbc_lnk_multi( divu_i (:,:), 'T', 1., delta_i(:,:), 'T', 1.,  shear_i(:,:), 'T', 1. ) 
    669  
    670       ! * Store the stress tensor for the next time step 
     581      CALL lbc_lnk_multi( delta_i, 'T', 1.,  shear_i, 'T', 1. ) 
     582       
     583      ! --- Store the stress tensor for the next time step --- ! 
    671584      stress1_i (:,:) = zs1 (:,:) 
    672585      stress2_i (:,:) = zs2 (:,:) 
    673586      stress12_i(:,:) = zs12(:,:) 
    674  
    675       ! 
     587      ! 
     588 
    676589      !------------------------------------------------------------------------------! 
    677590      ! 6) Control prints of residual and charge ellipse 
    … …  
    695608            WRITE(charout,FMT="('lim_rhg  :', I4, I6, I1, I1, A10)") 1000, numit, 0, 0, ' ch. ell. ' 
    696609            CALL prt_ctl_info(charout) 
    697             DO jj = k_j1+1, k_jpj-1 
     610            DO jj = 2, jpjm1 
    698611               DO ji = 2, jpim1 
    699612                  IF (zpresh(ji,jj) > 1.0) THEN 
    700                      sigma1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) )  
    701                      sigma2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 
     613                     zsig1 = ( zs1(ji,jj) + SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*zpresh(ji,jj) )  
     614                     zsig2 = ( zs1(ji,jj) - SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*zpresh(ji,jj) ) 
    702615                     WRITE(charout,FMT="('lim_rhg  :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)") 
    703616                     CALL prt_ctl_info(charout) 
    … …  
    710623      ENDIF 
    711624      ! 
    712       CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 
    713       CALL wrk_dealloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask               ) 
    714       CALL wrk_dealloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  ) 
    715       CALL wrk_dealloc( jpi,jpj, zs1   , zs2   , zs12   , zresr , zpice                 ) 
     625      CALL wrk_dealloc( jpi,jpj, zpresh, z1_e1t0, z1_e2t0, zp_delt ) 
     626      CALL wrk_dealloc( jpi,jpj, za1, za2, zmass1, zmass2, zcor1, zcor2 ) 
     627      CALL wrk_dealloc( jpi,jpj, zspgu, zspgv, v_oce1, u_oce2, v_ice1, u_ice2, zf1, zf2 ) 
     628      CALL wrk_dealloc( jpi,jpj, zds, zdt, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 
     629      CALL wrk_dealloc( jpi,jpj, zswitch1, zswitch2 ) 
    716630 
    717631   END SUBROUTINE lim_rhg 
    … …  
    722636   !!---------------------------------------------------------------------- 
    723637CONTAINS 
    724    SUBROUTINE lim_rhg( k1 , k2 )         ! Dummy routine 
    725       WRITE(*,*) 'lim_rhg: You should not have seen this print! error?', k1, k2 
     638   SUBROUTINE lim_rhg         ! Dummy routine 
     639      WRITE(*,*) 'lim_rhg: You should not have seen this print! error?' 
    726640   END SUBROUTINE lim_rhg 
    727641#endif 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r6515 r6746  
    2626   USE sbc_ice        ! Surface boundary condition: ice fields 
    2727   USE thd_ice        ! LIM thermodynamic sea-ice variables 
    28    USE dom_ice        ! LIM sea-ice domain 
    2928   USE limthd_dif     ! LIM: thermodynamics, vertical diffusion 
    3029   USE limthd_dh      ! LIM: thermodynamics, ice and snow thickness variation 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90

    r6316 r6746  
    2121   USE sbc_ice        ! Surface boundary condition: ice fields 
    2222   USE thd_ice        ! LIM thermodynamics 
    23    USE dom_ice        ! LIM domain 
    2423   USE ice            ! LIM variables 
    2524   USE limtab         ! LIM 2D <==> 1D 
    … …  
    7170      !!               update ht_s_1d, ht_i_1d and tbif_1d(:,:)       
    7271      !!------------------------------------------------------------------------ 
    73       INTEGER ::   ji,jj,jk,jl      ! dummy loop indices 
    74       INTEGER ::   nbpac            ! local integers  
    75       INTEGER ::   ii, ij, iter     !   -       - 
    76       REAL(wp)  ::   ztmelts, zdv, zfrazb, zweight, zde                         ! local scalars 
     72      INTEGER  ::   ji,jj,jk,jl      ! dummy loop indices 
     73      INTEGER  ::   nbpac            ! local integers  
     74      INTEGER  ::   ii, ij, iter     !   -       - 
     75      REAL(wp) ::   ztmelts, zdv, zfrazb, zweight, zde                          ! local scalars 
    7776      REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf                     !   -      - 
    7877      REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      - 
    … …  
    154153 
    155154      ! Default new ice thickness 
    156       WHERE( qlead(:,:) < 0._wp ) ; hicol = rn_hnewice 
    157       ELSEWHERE                   ; hicol = 0._wp 
     155      WHERE( qlead(:,:) < 0._wp ) ; hicol(:,:) = rn_hnewice 
     156      ELSEWHERE                   ; hicol(:,:) = 0._wp 
    158157      END WHERE 
    159158 
    … …  
    170169         zgamafr = 0.03 
    171170 
    172          DO jj = 2, jpj 
    173             DO ji = 2, jpi 
    174                IF ( qlead(ji,jj) < 0._wp ) THEN 
     171         DO jj = 2, jpjm1 
     172            DO ji = 2, jpim1 
     173               IF ( qlead(ji,jj) < 0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN ! activated if cooling and no landfast 
    175174                  !------------- 
    176175                  ! Wind stress 
    … …  
    195194                  !------------------- 
    196195                  ! C-grid ice velocity 
    197                   rswitch = MAX(  0._wp, SIGN( 1._wp , at_i(ji,jj) )  ) 
    198                   zvgx    = rswitch * ( u_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)  + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp 
    199                   zvgy    = rswitch * ( v_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)  + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp 
     196                  zvgx    = ( u_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)  + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp 
     197                  zvgy    = ( v_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)  + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp 
    200198 
    201199                  !----------------------------------- 
    … …  
    203201                  !----------------------------------- 
    204202                  ! absolute relative velocity 
    205                   zvrel2 = MAX(  ( zvfrx - zvgx ) * ( zvfrx - zvgx )   & 
    206                      &         + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15 * 0.15 ) 
     203                  rswitch      = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 
     204                  zvrel2       = MAX(  ( zvfrx - zvgx ) * ( zvfrx - zvgx )   & 
     205                     &               + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15 * 0.15 ) * rswitch 
    207206                  zvrel(ji,jj) = SQRT( zvrel2 ) 
    208207 
    … …  
    219218                     zfp = ( hicol(ji,jj) - zhicrit ) * ( 3.0 * hicol(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2 
    220219 
    221                      hicol(ji,jj) = hicol(ji,jj) - zf/zfp 
     220                     hicol(ji,jj) = hicol(ji,jj) - zf / MAX( zfp, epsi20 ) 
    222221                     iter = iter + 1 
    223222                  END DO 
    … …  
    228227         END DO  
    229228         !  
    230          CALL lbc_lnk( zvrel(:,:), 'T', 1. ) 
    231          CALL lbc_lnk( hicol(:,:), 'T', 1. ) 
     229         CALL lbc_lnk( zvrel, 'T', 1. ) 
     230         CALL lbc_lnk( hicol, 'T', 1. ) 
    232231 
    233232      ENDIF ! End of computation of frazil ice collection thickness 
    … …  
    240239      ! Select points for new ice formation 
    241240      !------------------------------------- 
    242       ! This occurs if open water energy budget is negative 
     241      ! This occurs if open water energy budget is negative (cooling) and there is no landfast ice 
    243242      nbpac = 0 
    244243      npac(:) = 0 
    … …  
    246245      DO jj = 1, jpj 
    247246         DO ji = 1, jpi 
    248             IF ( qlead(ji,jj)  <  0._wp ) THEN 
     247            IF ( qlead(ji,jj)  <  0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN 
    249248               nbpac = nbpac + 1 
    250249               npac( nbpac ) = (jj - 1) * jpi + ji 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90

    r6584 r6746  
    1717   USE dom_oce        ! ocean domain 
    1818   USE sbc_oce        ! ocean surface boundary condition 
    19    USE dom_ice        ! ice domain 
    2019   USE ice            ! ice variables 
    2120   USE limadv         ! ice advection 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90

    r6515 r6746  
    1515   USE sbc_oce         ! Surface boundary condition: ocean fields 
    1616   USE sbc_ice         ! Surface boundary condition: ice fields 
    17    USE dom_ice 
    1817   USE dom_oce 
    1918   USE phycst          ! physical constants 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90

    r6515 r6746  
    1515   USE sbc_oce         ! Surface boundary condition: ocean fields 
    1616   USE sbc_ice         ! Surface boundary condition: ice fields 
    17    USE dom_ice 
    1817   USE dom_oce 
    1918   USE phycst          ! physical constants 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90

    r6584 r6746  
    4141   USE ice            ! ice variables 
    4242   USE thd_ice        ! ice variables (thermodynamics) 
    43    USE dom_ice        ! ice domain 
    4443   USE in_out_manager ! I/O manager 
    4544   USE lib_mpp        ! MPP library 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90

    r6515 r6746  
    1717   USE sbc_oce         ! Surface boundary condition: ocean fields 
    1818   USE sbc_ice         ! Surface boundary condition: ice fields 
    19    USE dom_ice 
    2019   USE ice 
    2120   USE limvar 
    … …  
    4039   !!---------------------------------------------------------------------- 
    4140CONTAINS 
    42  
    43 #if defined key_dimgout 
    44 # include "limwri_dimg.h90" 
    45 #else 
    4641 
    4742   SUBROUTINE lim_wri( kindic ) 
    … …  
    5954      INTEGER  ::  ji, jj, jk, jl  ! dummy loop indices 
    6055      REAL(wp) ::  z1_365 
    61       REAL(wp) ::  ztmp 
     56      REAL(wp) ::  z2da, z2db, ztmp 
    6257      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zswi2 
    63       REAL(wp), POINTER, DIMENSION(:,:)   ::  z2d, z2da, z2db, zswi    ! 2D workspace 
     58      REAL(wp), POINTER, DIMENSION(:,:)   ::  z2d, zswi    ! 2D workspace 
    6459      !!------------------------------------------------------------------- 
    6560 
    … …  
    6762 
    6863      CALL wrk_alloc( jpi, jpj, jpl, zswi2 ) 
    69       CALL wrk_alloc( jpi, jpj     , z2d, z2da, z2db, zswi ) 
     64      CALL wrk_alloc( jpi, jpj     , z2d, zswi ) 
    7065 
    7166      !----------------------------- 
    … …  
    9590         DO jj = 2 , jpjm1 
    9691            DO ji = 2 , jpim1 
    97                z2da(ji,jj)  = (  u_ice(ji,jj) * umask(ji,jj,1) + u_ice(ji-1,jj) * umask(ji-1,jj,1) ) * 0.5_wp 
    98                z2db(ji,jj)  = (  v_ice(ji,jj) * vmask(ji,jj,1) + v_ice(ji,jj-1) * vmask(ji,jj-1,1) ) * 0.5_wp 
     92               z2da  = ( u_ice(ji,jj) * umask(ji,jj,1) + u_ice(ji-1,jj) * umask(ji-1,jj,1) ) * 0.5_wp 
     93               z2db  = ( v_ice(ji,jj) * vmask(ji,jj,1) + v_ice(ji,jj-1) * vmask(ji,jj-1,1) ) * 0.5_wp 
     94               z2d(ji,jj) = SQRT( z2da * z2da + z2db * z2db ) 
    9995           END DO 
    10096         END DO 
    101          CALL lbc_lnk( z2da, 'T', -1. ) 
    102          CALL lbc_lnk( z2db, 'T', -1. ) 
    103          CALL iom_put( "uice_ipa"     , z2da             )       ! ice velocity u component 
    104          CALL iom_put( "vice_ipa"     , z2db             )       ! ice velocity v component 
    105          DO jj = 1, jpj                                  
    106             DO ji = 1, jpi 
    107                z2d(ji,jj)  = SQRT( z2da(ji,jj) * z2da(ji,jj) + z2db(ji,jj) * z2db(ji,jj) )  
    108             END DO 
    109          END DO 
    110          CALL iom_put( "icevel"       , z2d * zswi       )       ! ice velocity module 
     97         CALL lbc_lnk( z2d, 'T', 1. ) 
     98         CALL iom_put( "uice_ipa"     , u_ice * zswi     )       ! ice velocity u component 
     99         CALL iom_put( "vice_ipa"     , v_ice * zswi     )       ! ice velocity v component 
     100         CALL iom_put( "icevel"       , z2d   * zswi     )       ! ice velocity module 
    111101      ENDIF 
    112102      ! 
    … …  
    130120      CALL iom_put( "micesalt"    , smt_i   * zswi      )        ! mean ice salinity 
    131121 
    132       CALL iom_put( "icestr"      , strength * 0.001 * zswi )    ! ice strength 
     122      CALL iom_put( "icestr"      , strength * zswi )    ! ice strength 
    133123      CALL iom_put( "idive"       , divu_i * 1.0e8   * zswi )    ! divergence 
    134124      CALL iom_put( "ishear"      , shear_i * 1.0e8  * zswi )    ! shear 
    … …  
    163153      CALL iom_put( "vfxlam"     , wfx_lam * ztmp       )        ! lateral melt  
    164154      CALL iom_put( "vfxice"     , wfx_ice * ztmp       )        ! total ice growth/melt  
     155 
     156      IF ( iom_use( "vfxthin" ) ) THEN   ! ice production for open water + thin ice (<20cm) => comparable to observations   
     157         WHERE( htm_i(:,:) < 0.2 .AND. htm_i(:,:) > 0. ) ; z2d = wfx_bog 
     158         ELSEWHERE                                       ; z2d = 0._wp 
     159         END WHERE 
     160         CALL iom_put( "vfxthin", ( wfx_opw + z2d ) * ztmp ) 
     161      ENDIF 
     162 
     163      ztmp = rday / rhosn 
     164      CALL iom_put( "vfxspr"     , wfx_spr * ztmp       )        ! precip (snow) 
    165165      CALL iom_put( "vfxsnw"     , wfx_snw * ztmp       )        ! total snw growth/melt  
    166       CALL iom_put( "vfxsub"     , wfx_sub * ztmp       )        ! sublimation (snow)  
    167       CALL iom_put( "vfxspr"     , wfx_spr * ztmp       )        ! precip (snow) 
     166      CALL iom_put( "vfxsub"     , wfx_sub * ztmp       )        ! sublimation (snow/ice)  
    168167       
    169168      CALL iom_put( "afxtot"     , afx_tot * rday       )        ! concentration tendency (total) 
    … …  
    190189      CALL iom_put ('hfxspr'     , hfx_spr(:,:)         )   ! Heat content of snow precip  
    191190 
    192  
    193       IF ( iom_use( "vfxthin" ) ) THEN   ! ice production for open water + thin ice (<20cm) => comparable to observations   
    194          WHERE( htm_i(:,:) < 0.2 .AND. htm_i(:,:) > 0. ) ; z2d = wfx_bog 
    195          ELSEWHERE                                       ; z2d = 0._wp 
    196          END WHERE 
    197          CALL iom_put( "vfxthin", ( wfx_opw + z2d ) * ztmp ) 
    198       ENDIF 
    199191       
    200192      !-------------------------------- 
    … …  
    223215       
    224216      CALL wrk_dealloc( jpi, jpj, jpl, zswi2 ) 
    225       CALL wrk_dealloc( jpi, jpj     , z2d, zswi, z2da, z2db ) 
     217      CALL wrk_dealloc( jpi, jpj     , z2d, zswi ) 
    226218 
    227219      IF( nn_timing == 1 )  CALL timing_stop('limwri') 
    228220       
    229221   END SUBROUTINE lim_wri 
    230 #endif 
    231222 
    232223  

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/NST_SRC/agrif_lim3_interp.F90

    r6584 r6746  
    2323   USE sbc_oce 
    2424   USE ice 
    25    USE dom_ice 
    2625   USE agrif_ice 
    2726 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/NST_SRC/agrif_lim3_update.F90

    r6584 r6746  
    2626   USE agrif_oce 
    2727   USE ice 
    28    USE dom_ice 
    2928   USE agrif_ice  
    3029 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90

    r6204 r6746  
    2727#elif defined key_lim3 
    2828   USE ice             ! LIM_3 ice variables 
    29    USE dom_ice         ! sea-ice domain 
    3029   USE limvar 
    3130#endif  

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90

    r5253 r6746  
    211211      REAL(wp) ::   zztmp   
    212212      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zsaldta   ! Jan/Dec levitus salinity 
    213       ! reading initial file 
    214       LOGICAL  ::   ln_tsd_init      !: T & S data flag 
    215       LOGICAL  ::   ln_tsd_tradmp    !: internal damping toward input data flag 
    216       CHARACTER(len=100)            ::   cn_dir 
    217       TYPE(FLD_N)                   ::  sn_tem,sn_sal 
    218       INTEGER  ::   ios=0 
    219  
    220       NAMELIST/namtsd/ ln_tsd_init,ln_tsd_tradmp,cn_dir,sn_tem,sn_sal 
    221       ! 
    222  
    223       REWIND( numnam_ref )              ! Namelist namtsd in reference namelist : 
    224       READ  ( numnam_ref, namtsd, IOSTAT = ios, ERR = 901) 
    225 901   IF( ios /= 0 ) CALL ctl_nam ( ios , ' namtsd in reference namelist for dia_ar5', lwp ) 
    226       REWIND( numnam_cfg )              ! Namelist namtsd in configuration namelist : Parameters of the run 
    227       READ  ( numnam_cfg, namtsd, IOSTAT = ios, ERR = 902 ) 
    228 902   IF( ios /= 0 ) CALL ctl_nam ( ios , ' namtsd in configuration namelist for dia_ar5', lwp ) 
    229       IF(lwm) WRITE ( numond, namtsd ) 
    230213      ! 
    231214      !!---------------------------------------------------------------------- 
    … …  
    233216      IF( nn_timing == 1 )   CALL timing_start('dia_ar5_init') 
    234217      ! 
    235       CALL wrk_alloc( jpi , jpj , jpk, jpts, zsaldta ) 
     218      CALL wrk_alloc( jpi, jpj, jpk, 2, zsaldta ) 
    236219      !                                      ! allocate dia_ar5 arrays 
    237220      IF( dia_ar5_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ar5_init : unable to allocate arrays' ) 
    … …  
    249232      IF( lk_mpp )   CALL mpp_sum( vol0 ) 
    250233 
    251       CALL iom_open ( TRIM( cn_dir )//TRIM(sn_sal%clname), inum ) 
    252       CALL iom_get  ( inum, jpdom_data, TRIM(sn_sal%clvar), zsaldta(:,:,:,1), 1  ) 
    253       CALL iom_get  ( inum, jpdom_data, TRIM(sn_sal%clvar), zsaldta(:,:,:,2), 12 ) 
     234      CALL iom_open ( 'sali_ref_clim_monthly', inum ) 
     235      CALL iom_get  ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,1), 1  ) 
     236      CALL iom_get  ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,2), 12 ) 
    254237      CALL iom_close( inum ) 
     238 
    255239      sn0(:,:,:) = 0.5_wp * ( zsaldta(:,:,:,1) + zsaldta(:,:,:,2) )         
    256240      sn0(:,:,:) = sn0(:,:,:) * tmask(:,:,:) 
    … …  
    267251      ENDIF 
    268252      ! 
    269       CALL wrk_dealloc( jpi , jpj , jpk, jpts, zsaldta ) 
     253      CALL wrk_dealloc( jpi, jpj, jpk, 2, zsaldta ) 
    270254      ! 
    271255      IF( nn_timing == 1 )   CALL timing_stop('dia_ar5_init') 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/LDF/ldfeiv.F90

    r4990 r6746  
    157157         END DO 
    158158      ENDIF 
     159 
     160      ! ORCA R1: Take the minimum between aeiw  and aeiv0 
     161      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN 
     162         DO jj = 2, jpjm1 
     163            DO ji = fs_2, fs_jpim1   ! vector opt. 
     164               aeiw(ji,jj) = MIN( aeiw(ji,jj), aeiv0 ) 
     165            END DO 
     166         END DO 
     167      ENDIF 
     168 
    159169      CALL lbc_lnk( aeiw, 'W', 1. )      ! lateral boundary condition on aeiw  
    160170 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r6584 r6746  
    13351335      !!             ***  ROUTINE sbc_cpl_ice_flx  *** 
    13361336      !! 
    1337       !! ** Purpose :   provide the heat and freshwater fluxes of the  
    1338       !!              ocean-ice system. 
     1337      !! ** Purpose :   provide the heat and freshwater fluxes of the ocean-ice system 
    13391338      !! 
    13401339      !! ** Method  :   transform the fields received from the atmosphere into 
    13411340      !!             surface heat and fresh water boundary condition for the  
    13421341      !!             ice-ocean system. The following fields are provided: 
    1343       !!              * total non solar, solar and freshwater fluxes (qns_tot,  
     1342      !!               * total non solar, solar and freshwater fluxes (qns_tot,  
    13441343      !!             qsr_tot and emp_tot) (total means weighted ice-ocean flux) 
    13451344      !!             NB: emp_tot include runoffs and calving. 
    1346       !!              * fluxes over ice (qns_ice, qsr_ice, emp_ice) where 
     1345      !!               * fluxes over ice (qns_ice, qsr_ice, emp_ice) where 
    13471346      !!             emp_ice = sublimation - solid precipitation as liquid 
    13481347      !!             precipitation are re-routed directly to the ocean and  
    1349       !!             runoffs and calving directly enter the ocean. 
    1350       !!              * solid precipitation (sprecip), used to add to qns_tot  
     1348      !!             calving directly enter the ocean (runoffs are read but included in trasbc.F90) 
     1349      !!               * solid precipitation (sprecip), used to add to qns_tot  
    13511350      !!             the heat lost associated to melting solid precipitation 
    13521351      !!             over the ocean fraction. 
    1353       !!       ===>> CAUTION here this changes the net heat flux received from 
    1354       !!             the atmosphere 
    1355       !! 
    1356       !!                  - the fluxes have been separated from the stress as 
    1357       !!                 (a) they are updated at each ice time step compare to 
    1358       !!                 an update at each coupled time step for the stress, and 
    1359       !!                 (b) the conservative computation of the fluxes over the 
    1360       !!                 sea-ice area requires the knowledge of the ice fraction 
    1361       !!                 after the ice advection and before the ice thermodynamics, 
    1362       !!                 so that the stress is updated before the ice dynamics 
    1363       !!                 while the fluxes are updated after it. 
     1352      !!               * heat content of rain, snow and evap can also be provided, 
     1353      !!             otherwise heat flux associated with these mass flux are 
     1354      !!             guessed (qemp_oce, qemp_ice) 
     1355      !! 
     1356      !!             - the fluxes have been separated from the stress as 
     1357      !!               (a) they are updated at each ice time step compare to 
     1358      !!               an update at each coupled time step for the stress, and 
     1359      !!               (b) the conservative computation of the fluxes over the 
     1360      !!               sea-ice area requires the knowledge of the ice fraction 
     1361      !!               after the ice advection and before the ice thermodynamics, 
     1362      !!               so that the stress is updated before the ice dynamics 
     1363      !!               while the fluxes are updated after it. 
     1364      !! 
     1365      !! ** Details 
     1366      !!             qns_tot = pfrld * qns_oce + ( 1 - pfrld ) * qns_ice   => provided 
     1367      !!                     + qemp_oce + qemp_ice                         => recalculated and added up to qns 
     1368      !! 
     1369      !!             qsr_tot = pfrld * qsr_oce + ( 1 - pfrld ) * qsr_ice   => provided 
     1370      !! 
     1371      !!             emp_tot = emp_oce + emp_ice                           => calving is provided and added to emp_tot (and emp_oce) 
     1372      !!                                                                      river runoff (rnf) is provided but not included here 
    13641373      !! 
    13651374      !! ** Action  :   update at each nf_ice time step: 
    13661375      !!                   qns_tot, qsr_tot  non-solar and solar total heat fluxes 
    13671376      !!                   qns_ice, qsr_ice  non-solar and solar heat fluxes over the ice 
    1368       !!                   emp_tot            total evaporation - precipitation(liquid and solid) (-runoff)(-calving) 
    1369       !!                   emp_ice            ice sublimation - solid precipitation over the ice 
    1370       !!                   dqns_ice           d(non-solar heat flux)/d(Temperature) over the ice 
    1371       !!                   sprecip             solid precipitation over the ocean   
     1377      !!                   emp_tot           total evaporation - precipitation(liquid and solid) (-calving) 
     1378      !!                   emp_ice           ice sublimation - solid precipitation over the ice 
     1379      !!                   dqns_ice          d(non-solar heat flux)/d(Temperature) over the ice 
     1380      !!                   sprecip           solid precipitation over the ocean   
    13721381      !!---------------------------------------------------------------------- 
    13731382      REAL(wp), INTENT(in   ), DIMENSION(:,:)   ::   p_frld     ! lead fraction                [0 to 1] 
    … …  
    13791388      INTEGER ::   jl         ! dummy loop index 
    13801389      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zcptn, ztmp, zicefr, zmsk, zsnw 
    1381       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap, zevap_ice, zdevap_ice 
     1390      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice 
    13821391      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice 
    13831392      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice 
    … …  
    13871396      ! 
    13881397      CALL wrk_alloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zsnw ) 
    1389       CALL wrk_alloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap, zevap_ice, zdevap_ice ) 
     1398      CALL wrk_alloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice ) 
    13901399      CALL wrk_alloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice ) 
    13911400      CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice ) 
    … …  
    13961405      ! 
    13971406      !                                                      ! ========================= ! 
    1398       !                                                      !    freshwater budget      !   (emp) 
     1407      !                                                      !    freshwater budget      !   (emp_tot) 
    13991408      !                                                      ! ========================= ! 
    14001409      ! 
    1401       !                                                           ! total Precipitation - total Evaporation (emp_tot) 
    1402       !                                                           ! solid precipitation - sublimation       (emp_ice) 
    1403       !                                                           ! solid Precipitation                     (sprecip) 
    1404       !                                                           ! liquid + solid Precipitation            (tprecip) 
     1410      !                                                           ! solid Precipitation                                (sprecip) 
     1411      !                                                           ! liquid + solid Precipitation                       (tprecip) 
     1412      !                                                           ! total Evaporation - total Precipitation            (emp_tot) 
     1413      !                                                           ! sublimation - solid precipitation (cell average)   (emp_ice) 
    14051414      SELECT CASE( TRIM( sn_rcv_emp%cldes ) ) 
    1406       CASE( 'conservative'  )   ! received fields: jpr_rain, jpr_snow, jpr_ievp, jpr_tevp 
    1407          zsprecip(:,:) = frcv(jpr_snow)%z3(:,:,1)                  ! May need to ensure positive here 
    1408          ztprecip(:,:) = frcv(jpr_rain)%z3(:,:,1) + zsprecip(:,:)  ! May need to ensure positive here 
    1409          zemp_tot(:,:) = frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:) 
    1410          zemp_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1) 
    1411             CALL iom_put( 'rain'         , frcv(jpr_rain)%z3(:,:,1)              )   ! liquid precipitation  
     1415      CASE( 'conservative' )   ! received fields: jpr_rain, jpr_snow, jpr_ievp, jpr_tevp 
     1416         zsprecip(:,:) =   frcv(jpr_snow)%z3(:,:,1)                  ! May need to ensure positive here 
     1417         ztprecip(:,:) =   frcv(jpr_rain)%z3(:,:,1) + zsprecip(:,:)  ! May need to ensure positive here 
     1418         zemp_tot(:,:) =   frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:) 
     1419         zemp_ice(:,:) = ( frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1) ) * zicefr(:,:) 
     1420               CALL iom_put( 'rain'         ,   frcv(jpr_rain)%z3(:,:,1)                                                         )  ! liquid precipitation  
    14121421         IF( iom_use('hflx_rain_cea') )   & 
    1413             CALL iom_put( 'hflx_rain_cea', frcv(jpr_rain)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from liq. precip.  
    1414          IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') )   & 
    1415             ztmp(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) 
     1422            &  CALL iom_put( 'hflx_rain_cea',   frcv(jpr_rain)%z3(:,:,1) * zcptn(:,:)                                            )  ! heat flux from liq. precip.  
    14161423         IF( iom_use('evap_ao_cea'  ) )   & 
    1417             CALL iom_put( 'evap_ao_cea'  , ztmp                   )   ! ice-free oce evap (cell average) 
     1424            &  CALL iom_put( 'evap_ao_cea'  ,   frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:)                )  ! ice-free oce evap (cell average) 
    14181425         IF( iom_use('hflx_evap_cea') )   & 
    1419             CALL iom_put( 'hflx_evap_cea', ztmp(:,:) * zcptn(:,:) )   ! heat flux from from evap (cell average) 
    1420       CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp 
     1426            &  CALL iom_put( 'hflx_evap_cea', ( frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) ) * zcptn(:,:) )  ! heat flux from from evap (cell average) 
     1427      CASE( 'oce and ice' )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp 
    14211428         zemp_tot(:,:) = p_frld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + zicefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1) 
    1422          zemp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1) 
     1429         zemp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1) * zicefr(:,:) 
    14231430         zsprecip(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_semp)%z3(:,:,1) 
    14241431         ztprecip(:,:) = frcv(jpr_semp)%z3(:,:,1) - frcv(jpr_sbpr)%z3(:,:,1) + zsprecip(:,:) 
    … …  
    14261433 
    14271434#if defined key_lim3 
    1428       ! zsnw = snow percentage over ice after wind blowing 
    1429       zsnw(:,:) = 0._wp 
    1430       CALL lim_thd_snwblow( p_frld, zsnw ) 
     1435      ! zsnw = snow fraction over ice after wind blowing 
     1436      zsnw(:,:) = 0._wp  ;  CALL lim_thd_snwblow( p_frld, zsnw ) 
    14311437       
    1432       ! --- evaporation (kg/m2/s) --- ! 
     1438      ! --- evaporation minus precipitation corrected (because of wind blowing on snow) --- ! 
     1439      zemp_ice(:,:) = zemp_ice(:,:) + zsprecip(:,:) * ( zicefr(:,:) - zsnw(:,:) )  ! emp_ice = A * sublimation - zsnw * sprecip 
     1440      zemp_oce(:,:) = zemp_tot(:,:) - zemp_ice(:,:)                                ! emp_oce = emp_tot - emp_ice 
     1441 
     1442      ! --- evaporation over ocean (used later for qemp) --- ! 
     1443      zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) 
     1444 
     1445      ! --- evaporation over ice (kg/m2/s) --- ! 
    14331446      zevap_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1) 
    14341447      ! since the sensitivity of evap to temperature (devap/dT) is not prescribed by the atmosphere, we set it to 0 
    … …  
    14361449      zdevap_ice(:,:) = 0._wp 
    14371450       
    1438       ! --- evaporation minus precipitation corrected for the effect of wind blowing on snow --- ! 
    1439       zemp_oce(:,:) = zemp_tot(:,:) - zemp_ice(:,:) - zsprecip * (1._wp - zsnw) 
    1440       zemp_ice(:,:) = zemp_ice(:,:) + zsprecip * (1._wp - zsnw)           
    1441  
    1442       ! Sublimation over sea-ice (cell average) 
    1443       IF( iom_use('subl_ai_cea') )  CALL iom_put( 'subl_ai_cea', zevap_ice(:,:) * zicefr(:,:) ) 
    1444       ! runoffs and calving (put in emp_tot) 
     1451      ! --- runoffs (included in emp later on) --- ! 
    14451452      IF( srcv(jpr_rnf)%laction )   rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1) 
     1453 
     1454      ! --- calving (put in emp_tot and emp_oce) --- ! 
    14461455      IF( srcv(jpr_cal)%laction ) THEN  
    14471456         zemp_tot(:,:) = zemp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1) 
     1457         zemp_oce(:,:) = zemp_oce(:,:) - frcv(jpr_cal)%z3(:,:,1) 
    14481458         CALL iom_put( 'calving_cea', frcv(jpr_cal)%z3(:,:,1) ) 
    14491459      ENDIF 
    … …  
    14711481      ENDIF 
    14721482 
    1473                                      CALL iom_put( 'snowpre'    , sprecip                         )  ! Snow 
    1474       IF( iom_use('snow_ao_cea') )   CALL iom_put( 'snow_ao_cea', sprecip(:,:) * ( 1._wp - zsnw ) )  ! Snow over ice-free ocean  (cell average) 
    1475       IF( iom_use('snow_ai_cea') )   CALL iom_put( 'snow_ai_cea', sprecip(:,:) *           zsnw   )  ! Snow over sea-ice         (cell average)     
     1483      IF( iom_use('subl_ai_cea') )   CALL iom_put( 'subl_ai_cea', zevap_ice(:,:) * zicefr(:,:)         )  ! Sublimation over sea-ice (cell average) 
     1484                                     CALL iom_put( 'snowpre'    , sprecip(:,:)                         )  ! Snow 
     1485      IF( iom_use('snow_ao_cea') )   CALL iom_put( 'snow_ao_cea', sprecip(:,:) * ( 1._wp - zsnw(:,:) ) )  ! Snow over ice-free ocean  (cell average) 
     1486      IF( iom_use('snow_ai_cea') )   CALL iom_put( 'snow_ai_cea', sprecip(:,:) *           zsnw(:,:)   )  ! Snow over sea-ice         (cell average) 
    14761487#else 
    1477       ! Sublimation over sea-ice (cell average) 
    1478       IF( iom_use('subl_ai_cea') )  CALL iom_put( 'subl_ai_cea', frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) ) 
    14791488      ! runoffs and calving (put in emp_tot) 
    14801489      IF( srcv(jpr_rnf)%laction )   rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1) 
    … …  
    14961505      ENDIF 
    14971506 
    1498          CALL iom_put( 'snowpre'    , sprecip                                )   ! Snow 
    1499       IF( iom_use('snow_ao_cea') )   & 
    1500          CALL iom_put( 'snow_ao_cea', sprecip(:,:) * p_frld(:,:)             )   ! Snow        over ice-free ocean  (cell average) 
    1501       IF( iom_use('snow_ai_cea') )   & 
    1502          CALL iom_put( 'snow_ai_cea', sprecip(:,:) * zicefr(:,:)             )   ! Snow        over sea-ice         (cell average) 
     1507      IF( iom_use('subl_ai_cea') )  CALL iom_put( 'subl_ai_cea', frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) )  ! Sublimation over sea-ice (cell average) 
     1508                                    CALL iom_put( 'snowpre'    , sprecip(:,:)               )   ! Snow 
     1509      IF( iom_use('snow_ao_cea') )  CALL iom_put( 'snow_ao_cea', sprecip(:,:) * p_frld(:,:) )   ! Snow over ice-free ocean  (cell average) 
     1510      IF( iom_use('snow_ai_cea') )  CALL iom_put( 'snow_ai_cea', sprecip(:,:) * zicefr(:,:) )   ! Snow over sea-ice         (cell average) 
    15031511#endif 
    15041512 
    … …  
    15061514      SELECT CASE( TRIM( sn_rcv_qns%cldes ) )                !   non solar heat fluxes   !   (qns) 
    15071515      !                                                      ! ========================= ! 
    1508       CASE( 'oce only' )                                     ! the required field is directly provided 
    1509          zqns_tot(:,:  ) = frcv(jpr_qnsoce)%z3(:,:,1) 
    1510       CASE( 'conservative' )                                      ! the required fields are directly provided 
    1511          zqns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1) 
     1516      CASE( 'oce only' )         ! the required field is directly provided 
     1517         zqns_tot(:,:) = frcv(jpr_qnsoce)%z3(:,:,1) 
     1518      CASE( 'conservative' )     ! the required fields are directly provided 
     1519         zqns_tot(:,:) = frcv(jpr_qnsmix)%z3(:,:,1) 
    15121520         IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN 
    15131521            zqns_ice(:,:,1:jpl) = frcv(jpr_qnsice)%z3(:,:,1:jpl) 
    15141522         ELSE 
    1515             ! Set all category values equal for the moment 
    15161523            DO jl=1,jpl 
    1517                zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1) 
     1524               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1) ! Set all category values equal 
    15181525            ENDDO 
    15191526         ENDIF 
    1520       CASE( 'oce and ice' )       ! the total flux is computed from ocean and ice fluxes 
    1521          zqns_tot(:,:  ) =  p_frld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1) 
     1527      CASE( 'oce and ice' )      ! the total flux is computed from ocean and ice fluxes 
     1528         zqns_tot(:,:) =  p_frld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1) 
    15221529         IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN 
    15231530            DO jl=1,jpl 
    … …  
    15261533            ENDDO 
    15271534         ELSE 
    1528             qns_tot(:,:   ) = qns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1) 
     1535            qns_tot(:,:) = qns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1) 
    15291536            DO jl=1,jpl 
    15301537               zqns_tot(:,:   ) = zqns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1) 
    … …  
    15321539            ENDDO 
    15331540         ENDIF 
    1534       CASE( 'mixed oce-ice' )     ! the ice flux is cumputed from the total flux, the SST and ice informations 
     1541      CASE( 'mixed oce-ice' )    ! the ice flux is cumputed from the total flux, the SST and ice informations 
    15351542! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED ** 
    15361543         zqns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1) 
    15371544         zqns_ice(:,:,1) = frcv(jpr_qnsmix)%z3(:,:,1)    & 
    15381545            &            + frcv(jpr_dqnsdt)%z3(:,:,1) * ( pist(:,:,1) - ( (rt0 + psst(:,:  ) ) * p_frld(:,:)   & 
    1539             &                                                   +          pist(:,:,1)  * zicefr(:,:) ) ) 
     1546            &                                           + pist(:,:,1) * zicefr(:,:) ) ) 
    15401547      END SELECT 
    15411548!!gm 
    … …  
    15471554!! similar job should be done for snow and precipitation temperature 
    15481555      !                                      
    1549       IF( srcv(jpr_cal)%laction ) THEN                            ! Iceberg melting  
    1550          ztmp(:,:) = frcv(jpr_cal)%z3(:,:,1) * lfus               ! add the latent heat of iceberg melting  
    1551          zqns_tot(:,:) = zqns_tot(:,:) - ztmp(:,:) 
    1552          IF( iom_use('hflx_cal_cea') )   & 
    1553             CALL iom_put( 'hflx_cal_cea', ztmp + frcv(jpr_cal)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from calving 
    1554       ENDIF 
    1555  
    1556       ztmp(:,:) = p_frld(:,:) * zsprecip(:,:) * lfus 
    1557       IF( iom_use('hflx_snow_cea') )    CALL iom_put( 'hflx_snow_cea', ztmp + sprecip(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average) 
     1556      IF( srcv(jpr_cal)%laction ) THEN   ! Iceberg melting  
     1557         zqns_tot(:,:) = zqns_tot(:,:) - frcv(jpr_cal)%z3(:,:,1) * lfus  ! add the latent heat of iceberg melting 
     1558                                                                         ! we suppose it melts at 0deg, though it should be temp. of surrounding ocean 
     1559         IF( iom_use('hflx_cal_cea') )   CALL iom_put( 'hflx_cal_cea', - frcv(jpr_cal)%z3(:,:,1) * lfus )   ! heat flux from calving 
     1560      ENDIF 
    15581561 
    15591562#if defined key_lim3       
    1560       ! --- evaporation --- ! 
    1561       zevap(:,:) = zemp_tot(:,:) + ztprecip(:,:) ! evaporation over ocean 
    1562  
    15631563      ! --- non solar flux over ocean --- ! 
    15641564      !         note: p_frld cannot be = 0 since we limit the ice concentration to amax 
    … …  
    15671567 
    15681568      ! --- heat flux associated with emp (W/m2) --- ! 
    1569       zqemp_oce(:,:) = -      zevap(:,:)                   * p_frld(:,:)      *   zcptn(:,:)   &       ! evap 
     1569      zqemp_oce(:,:) = -  zevap_oce(:,:)                                      *   zcptn(:,:)   &       ! evap 
    15701570         &             + ( ztprecip(:,:) - zsprecip(:,:) )                    *   zcptn(:,:)   &       ! liquid precip 
    15711571         &             +   zsprecip(:,:)                   * ( 1._wp - zsnw ) * ( zcptn(:,:) - lfus )  ! solid precip over ocean + snow melting 
    … …  
    16061606         qemp_ice (:,:  ) = zqemp_ice (:,:  ) 
    16071607      ENDIF 
     1608 
     1609      !! clem: we should output qemp_oce and qemp_ice (at least) 
     1610      IF( iom_use('hflx_snow_cea') )   CALL iom_put( 'hflx_snow_cea', sprecip(:,:) * ( zcptn(:,:) - Lfus ) )   ! heat flux from snow (cell average) 
     1611      !! these diags are not outputed yet 
     1612!!      IF( iom_use('hflx_rain_cea') )   CALL iom_put( 'hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * zcptn(:,:) )   ! heat flux from rain (cell average) 
     1613!!      IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put( 'hflx_snow_ao_cea', sprecip(:,:) * ( zcptn(:,:) - Lfus ) * (1._wp - zsnw(:,:)) ) ! heat flux from snow (cell average) 
     1614!!      IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put( 'hflx_snow_ai_cea', sprecip(:,:) * ( zcptn(:,:) - Lfus ) * zsnw(:,:) ) ! heat flux from snow (cell average) 
     1615 
    16081616#else 
    16091617      ! clem: this formulation is certainly wrong... but better than it was... 
    … …  
    16111619         &          - ztmp(:,:)                           &            ! remove the latent heat flux of solid precip. melting 
    16121620         &          - (  zemp_tot(:,:)                    &            ! remove the heat content of mass flux (assumed to be at SST) 
    1613          &             - zemp_ice(:,:) * zicefr(:,:)  ) * zcptn(:,:)  
     1621         &             - zemp_ice(:,:) ) * zcptn(:,:)  
    16141622 
    16151623     IF( ln_mixcpl ) THEN 
    … …  
    17311739 
    17321740      CALL wrk_dealloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zsnw ) 
    1733       CALL wrk_dealloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap, zevap_ice, zdevap_ice ) 
     1741      CALL wrk_dealloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice ) 
    17341742      CALL wrk_dealloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice ) 
    17351743      CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice ) 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r6584 r6746  
    2424   USE ice             ! LIM-3: ice variables 
    2525   USE thd_ice         ! LIM-3: thermodynamical variables 
    26    USE dom_ice         ! LIM-3: ice domain 
    2726 
    2827   USE sbc_oce         ! Surface boundary condition: ocean fields 
    … …  
    4847   USE limvar          ! Ice variables switch 
    4948 
    50    USE limmsh          ! LIM mesh 
    5149   USE limistate       ! LIM initial state 
    5250   USE limthd_sal      ! LIM ice thermodynamics: salinity 
    … …  
    310308      !                                ! Allocate the ice arrays 
    311309      ierr =        ice_alloc        ()      ! ice variables 
    312       ierr = ierr + dom_ice_alloc    ()      ! domain 
    313310      ierr = ierr + sbc_ice_alloc    ()      ! surface forcing 
    314311      ierr = ierr + thd_ice_alloc    ()      ! thermodynamics 
    … …  
    331328      ! 
    332329      CALL lim_thd_sal_init            ! set ice salinity parameters 
    333       ! 
    334       CALL lim_msh                     ! ice mesh initialization 
    335330      ! 
    336331      IF( ln_limdyn )   CALL lim_itd_me_init             ! ice thickness distribution initialization for mecanical deformation 
    … …  
    661656      diag_heat(:,:) = 0._wp ;   diag_smvi(:,:) = 0._wp ; 
    662657      diag_vice(:,:) = 0._wp ;   diag_vsnw(:,:) = 0._wp ; 
     658 
     659      tau_icebfr(:,:) = 0._wp; ! landfast ice param only (clem: important to keep the init here) 
    663660       
    664661   END SUBROUTINE sbc_lim_diag0 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90

    r6204 r6746  
    173173            DO jj = 2, jpjm1 
    174174               DO ji = fs_2, fs_jpim1   ! vector opt. 
    175                   zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    176175                  ! total intermediate advective trends 
    177                   ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
    178                      &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
    179                      &             + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) 
     176                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
     177                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
     178                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) / e1e2t(ji,jj) 
    180179                  ! update and guess with monotonic sheme 
    181                   pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn)         + ztra  * tmask(ji,jj,jk) 
    182                   zwi(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk) 
     180                  pta(ji,jj,jk,jn) =                       pta(ji,jj,jk,jn) +         ztra   / fse3t_n(ji,jj,jk) * tmask(ji,jj,jk) 
     181                  zwi(ji,jj,jk)    = ( fse3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + z2dtt * ztra ) / fse3t_a(ji,jj,jk) * tmask(ji,jj,jk) 
    183182               END DO 
    184183            END DO 
    … …  
    410409            DO jj = 2, jpjm1 
    411410               DO ji = fs_2, fs_jpim1   ! vector opt. 
    412                   zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    413411                  ! total intermediate advective trends 
    414                   ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
    415                      &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
    416                      &             + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) 
     412                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
     413                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
     414                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) / e1e2t(ji,jj) 
    417415                  ! update and guess with monotonic sheme 
    418                   pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn)         + ztra 
    419                   zwi(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk) 
     416                  pta(ji,jj,jk,jn) =                       pta(ji,jj,jk,jn) +         ztra   / fse3t_n(ji,jj,jk) * tmask(ji,jj,jk) 
     417                  zwi(ji,jj,jk)    = ( fse3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + z2dtt * ztra ) / fse3t_a(ji,jj,jk) * tmask(ji,jj,jk) 
    420418               END DO 
    421419            END DO 
    … …  
    438436         ! -------------------------------------------------- 
    439437         ! antidiffusive flux on i and j 
    440  
    441  
    442          DO jk = 1, jpkm1 
    443  
     438         ! 
     439         DO jk = 1, jpkm1 
     440            ! 
    444441            DO jj = 1, jpjm1 
    445442               DO ji = 1, fs_jpim1   ! vector opt. 
    … …  
    572569   END SUBROUTINE tra_adv_tvd_zts 
    573570 
     571 
    574572   SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt ) 
    575573      !!--------------------------------------------------------------------- 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90

    r6471 r6746  
    158158         ELSE                                         ! No restart or restart not found: Euler forward time stepping 
    159159            zfact = 1._wp 
     160            sbc_tsc(:,:,:) = 0._wp 
    160161            sbc_tsc_b(:,:,:) = 0._wp 
    161162         ENDIF 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90

    r6204 r6746  
    7676      REAL(wp) ::   zchl 
    7777      REAL(wp) ::   zc0 , zc1 , zc2, zc3, z1_dep 
    78       REAL(wp), POINTER, DIMENSION(:,:  ) :: zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4, zqsr100 
     78      REAL(wp), POINTER, DIMENSION(:,:  ) :: zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 
     79      REAL(wp), POINTER, DIMENSION(:,:  ) :: zqsr100, zqsr_corr 
    7980      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpar, ze0, ze1, ze2, ze3 
    8081      !!--------------------------------------------------------------------- 
    … …  
    8384      ! 
    8485      ! Allocate temporary workspace 
    85       CALL wrk_alloc( jpi, jpj,      zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 ) 
     86      CALL wrk_alloc( jpi, jpj,      zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 ) 
     87      CALL wrk_alloc( jpi, jpj,      zqsr100, zqsr_corr ) 
    8688      CALL wrk_alloc( jpi, jpj, jpk, zpar, ze0, ze1, ze2, ze3 ) 
    8789 
    … …  
    112114      !                                        !  -------------------------------------- 
    113115      IF( l_trcdm2dc ) THEN                     !  diurnal cycle 
    114          ! 1% of qsr to compute euphotic layer 
     116         !                                       ! 1% of qsr to compute euphotic layer 
    115117         zqsr100(:,:) = 0.01 * qsr_mean(:,:)     !  daily mean qsr 
    116118         ! 
    117          CALL p4z_opt_par( kt, qsr_mean, ze1, ze2, ze3 )  
     119         zqsr_corr(:,:) = qsr_mean(:,:) / ( 1. - fr_i(:,:) + rtrn ) 
     120         ! 
     121         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 )  
    118122         ! 
    119123         DO jk = 1, nksrp       
    … …  
    123127         END DO 
    124128         ! 
    125          CALL p4z_opt_par( kt, qsr, ze1, ze2, ze3 )  
     129         zqsr_corr(:,:) = qsr(:,:) / ( 1. - fr_i(:,:) + rtrn ) 
     130         ! 
     131         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 )  
    126132         ! 
    127133         DO jk = 1, nksrp       
    … …  
    133139         zqsr100(:,:) = 0.01 * qsr(:,:) 
    134140         ! 
    135          CALL p4z_opt_par( kt, qsr, ze1, ze2, ze3 )  
     141         zqsr_corr(:,:) = qsr(:,:) / ( 1. - fr_i(:,:) + rtrn ) 
     142         ! 
     143         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 )  
    136144         ! 
    137145         DO jk = 1, nksrp       
    … …  
    226234      ENDIF 
    227235      ! 
    228       CALL wrk_dealloc( jpi, jpj,      zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 ) 
     236      CALL wrk_dealloc( jpi, jpj,      zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 ) 
     237      CALL wrk_dealloc( jpi, jpj,      zqsr100, zqsr_corr ) 
    229238      CALL wrk_dealloc( jpi, jpj, jpk, zpar,  ze0, ze1, ze2, ze3 ) 
    230239      ! 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90

    r6204 r6746  
    136136                  zval = MAX( 1., zstrn(ji,jj) ) 
    137137                  zval = 1.5 * zval / ( 12. + zval ) 
    138                   zprbio(ji,jj,jk) = prmax(ji,jj,jk) * zval 
     138                  zprbio(ji,jj,jk) = prmax(ji,jj,jk) * zval * ( 1. - fr_i(ji,jj) ) 
    139139                  zprdia(ji,jj,jk) = zprbio(ji,jj,jk) 
    140140               ENDIF 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/TRP/trcdmp.F90

    r6308 r6746  
    3535   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   restotr   ! restoring coeff. on tracers (s-1) 
    3636 
    37    INTEGER, PARAMETER           ::   npncts   = 5        ! number of closed sea 
     37   INTEGER, PARAMETER           ::   npncts   = 8        ! number of closed sea 
    3838   INTEGER, DIMENSION(npncts)   ::   nctsi1, nctsj1      ! south-west closed sea limits (i,j) 
    3939   INTEGER, DIMENSION(npncts)   ::   nctsi2, nctsj2      ! north-east closed sea limits (i,j) 
    … …  
    107107                
    108108               jl = n_trc_index(jn)  
    109                CALL trc_dta( kt, sf_trcdta(jl) )   ! read tracer data at nit000 
    110                ztrcdta(:,:,:) = sf_trcdta(jl)%fnow(:,:,:) * tmask(:,:,:) * rf_trfac(jl) 
     109               CALL trc_dta( kt, sf_trcdta(jl), rf_trfac(jl), ztrcdta )   ! read tracer data at nit000 
    111110 
    112111               SELECT CASE ( nn_zdmp_tr ) 
    … …  
    187186      INTEGER :: ji , jj, jk, jn, jl, jc                     ! dummy loop indicesa 
    188187      INTEGER :: isrow                                      ! local index 
     188      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrcdta       ! 3D  workspace 
    189189 
    190190      !!---------------------------------------------------------------------- 
    … …  
    207207            ! 
    208208                                                        ! Caspian Sea 
    209             nctsi1(1)   = 332  ; nctsj1(1)   = 243 - isrow 
    210             nctsi2(1)   = 344  ; nctsj2(1)   = 275 - isrow 
     209            nctsi1(1)   = 333  ; nctsj1(1)   = 243 - isrow 
     210            nctsi2(1)   = 342  ; nctsj2(1)   = 274 - isrow 
     211            !                                           ! Lake Superior 
     212            nctsi1(2)   = 198  ; nctsj1(2)   = 258 - isrow 
     213            nctsi2(2)   = 204  ; nctsj2(2)   = 262 - isrow 
     214            !                                           ! Lake Michigan 
     215            nctsi1(3)   = 201  ; nctsj1(3)   = 250 - isrow 
     216            nctsi2(3)   = 203  ; nctsj2(3)   = 256 - isrow 
     217            !                                           ! Lake Huron 
     218            nctsi1(4)   = 204  ; nctsj1(4)   = 252 - isrow 
     219            nctsi2(4)   = 209  ; nctsj2(4)   = 256 - isrow 
     220            !                                           ! Lake Erie 
     221            nctsi1(5)   = 206  ; nctsj1(5)   = 249 - isrow 
     222            nctsi2(5)   = 209  ; nctsj2(5)   = 251 - isrow 
     223            !                                           ! Lake Ontario 
     224            nctsi1(6)   = 210  ; nctsj1(6)   = 252 - isrow 
     225            nctsi2(6)   = 212  ; nctsj2(6)   = 252 - isrow 
     226            !                                           ! Victoria Lake 
     227            nctsi1(7)   = 321  ; nctsj1(7)   = 180 - isrow 
     228            nctsi2(7)   = 322  ; nctsj2(7)   = 189 - isrow 
     229            !                                           ! Baltic Sea 
     230            nctsi1(8)   = 297  ; nctsj1(8)   = 270 - isrow 
     231            nctsi2(8)   = 308  ; nctsj2(8)   = 293 - isrow 
    211232            !                                         
    212233            !                                           ! ======================= 
    … …  
    277298         IF(lwp)  WRITE(numout,*) 
    278299         ! 
     300         CALL wrk_alloc( jpi, jpj, jpk, ztrcdta )   ! Memory allocation 
     301         ! 
    279302         DO jn = 1, jptra 
    280303            IF( ln_trc_ini(jn) ) THEN      ! update passive tracers arrays with input data read from file 
    281304                jl = n_trc_index(jn) 
    282                 CALL trc_dta( kt, sf_trcdta(jl) )   ! read tracer data at nit000 
     305                CALL trc_dta( kt, sf_trcdta(jl), rf_trfac(jl), ztrcdta )   ! read tracer data at nit000 
    283306                DO jc = 1, npncts 
    284307                   DO jk = 1, jpkm1 
    285308                      DO jj = nctsj1(jc), nctsj2(jc) 
    286309                         DO ji = nctsi1(jc), nctsi2(jc) 
    287                             trn(ji,jj,jk,jn) = sf_trcdta(jl)%fnow(ji,jj,jk) * tmask(ji,jj,jk) * rf_trfac(jl) 
     310                            trn(ji,jj,jk,jn) = ztrcdta(ji,jj,jk) 
    288311                            trb(ji,jj,jk,jn) = trn(ji,jj,jk,jn) 
    289312                         ENDDO 
    … …  
    293316             ENDIF 
    294317          ENDDO 
    295           ! 
     318          CALL wrk_dealloc( jpi, jpj, jpk, ztrcdta ) 
    296319      ENDIF 
    297320      ! 
    … …  
    313336      IF( nn_timing == 1 )  CALL timing_start('trc_dmp_init') 
    314337      ! 
     338      !Allocate arrays 
     339      IF( trc_dmp_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'trc_dmp_init: unable to allocate arrays' ) 
    315340 
    316341      IF( lzoom )   nn_zdmp_tr = 0           ! restoring to climatology at closed north or south boundaries 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/trcdta.F90

    r6308 r6746  
    7777      ALLOCATE( n_trc_index(ntrc), slf_i(ntrc), STAT=ierr0 ) 
    7878      IF( ierr0 > 0 ) THEN 
    79          CALL ctl_stop( 'trc_nam: unable to allocate n_trc_index' )   ;   RETURN 
     79         CALL ctl_stop( 'trc_dta_init: unable to allocate n_trc_index' )   ;   RETURN 
    8080      ENDIF 
    8181      nb_trcdta      = 0 
    … …  
    9191      IF(lwp) THEN 
    9292         WRITE(numout,*) ' ' 
     93         WRITE(numout,*) 'trc_dta_init : Passive tracers Initial Conditions ' 
     94         WRITE(numout,*) '~~~~~~~~~~~~~~ ' 
    9395         WRITE(numout,*) ' number of passive tracers to be initialize by data :', ntra 
    9496         WRITE(numout,*) ' ' 
    … …  
    107109         DO jn = 1, ntrc 
    108110            IF( ln_trc_ini(jn) )  THEN    ! open input file only if ln_trc_ini(jn) is true 
    109                clndta = TRIM( sn_trcdta(jn)%clvar )  
    110                clntrc = TRIM( ctrcnm   (jn)       )  
     111               clndta = TRIM( sn_trcdta(jn)%clvar ) 
     112               if (jn > jptra) then 
     113                  clntrc='Dummy' ! By pass weird formats in ocean.output if ntrc > jptra 
     114               else 
     115                  clntrc = TRIM( ctrcnm   (jn)       ) 
     116               endif 
    111117               zfact  = rn_trfac(jn) 
    112                IF( clndta /=  clntrc ) THEN  
    113                   CALL ctl_warn( 'trc_dta_init: passive tracer data initialisation :  ',   & 
    114                   &              'the variable name in the data file : '//clndta//   &  
    115                   &              '  must be the same than the name of the passive tracer : '//clntrc//' ') 
     118               IF( clndta /=  clntrc ) THEN 
     119                  CALL ctl_warn( 'trc_dta_init: passive tracer data initialisation    ',   & 
     120                  &              'Input name of data file : '//TRIM(clndta)//   & 
     121                  &              ' differs from that of tracer : '//TRIM(clntrc)//' ') 
    116122               ENDIF 
    117                WRITE(numout,*) ' read an initial file for passive tracer number :', jn, ' name : ', clndta, &  
    118                &               ' multiplicative factor : ', zfact 
     123               WRITE(numout,'(a, i4,3a,e11.3)') ' Read IC file for tracer number :', & 
     124               &            jn, ', name : ', TRIM(clndta), ', Multiplicative Scaling factor : ', zfact 
    119125            ENDIF 
    120126         END DO 
    … …  
    124130         ALLOCATE( sf_trcdta(nb_trcdta), rf_trfac(nb_trcdta), STAT=ierr1 ) 
    125131         IF( ierr1 > 0 ) THEN 
    126             CALL ctl_stop( 'trc_dta_ini: unable to allocate  sf_trcdta structure' )   ;   RETURN 
     132            CALL ctl_stop( 'trc_dta_init: unable to allocate  sf_trcdta structure' )   ;   RETURN 
    127133         ENDIF 
    128134         ! 
    … …  
    135141               IF( sn_trcdta(jn)%ln_tint )  ALLOCATE( sf_trcdta(jl)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 ) 
    136142               IF( ierr2 + ierr3 > 0 ) THEN 
    137                  CALL ctl_stop( 'trc_dta : unable to allocate passive tracer data arrays' )   ;   RETURN 
     143                 CALL ctl_stop( 'trc_dta_init : unable to allocate passive tracer data arrays' )   ;   RETURN 
    138144               ENDIF 
    139145            ENDIF 
    … …  
    141147         ENDDO 
    142148         !                         ! fill sf_trcdta with slf_i and control print 
    143          CALL fld_fill( sf_trcdta, slf_i, cn_dir, 'trc_dta', 'Passive tracer data', 'namtrc' ) 
     149         CALL fld_fill( sf_trcdta, slf_i, cn_dir, 'trc_dta_init', 'Passive tracer data', 'namtrc' ) 
    144150         ! 
    145151      ENDIF 
    … …  
    151157 
    152158 
    153    SUBROUTINE trc_dta( kt, sf_dta ) 
     159   SUBROUTINE trc_dta( kt, sf_dta, ptrfac, ptrc) 
    154160      !!---------------------------------------------------------------------- 
    155161      !!                   ***  ROUTINE trc_dta  *** 
    … …  
    164170      !!---------------------------------------------------------------------- 
    165171      INTEGER                     , INTENT(in   ) ::   kt     ! ocean time-step 
    166       TYPE(FLD), DIMENSION(1)   , INTENT(inout) ::   sf_dta     ! array of information on the field to read 
     172      TYPE(FLD), DIMENSION(1)     , INTENT(inout) ::   sf_dta     ! array of information on the field to read 
     173      REAL(wp)                    , INTENT(in   ) ::   ptrfac  ! multiplication factor 
     174      REAL(wp), DIMENSION(jpi,jpj,jpk), OPTIONAL  , INTENT(out  ) ::   ptrc 
    167175      ! 
    168176      INTEGER ::   ji, jj, jk, jl, jkk, ik    ! dummy loop indices 
    169177      REAL(wp)::   zl, zi 
    170178      REAL(wp), DIMENSION(jpk) ::  ztp                ! 1D workspace 
     179      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrcdta   ! 3D  workspace 
    171180      CHARACTER(len=100) :: clndta 
    172181      !!---------------------------------------------------------------------- 
    … …  
    176185      IF( nb_trcdta > 0 ) THEN 
    177186         ! 
     187         CALL wrk_alloc( jpi, jpj, jpk, ztrcdta )    ! Memory allocation 
     188         ! 
    178189         CALL fld_read( kt, 1, sf_dta )      !==   read data at kt time step   ==! 
     190         ztrcdta(:,:,:) = sf_dta(1)%fnow(:,:,:) * tmask(:,:,:)    ! Mask 
    179191         ! 
    180192         IF( ln_sco ) THEN                   !==   s- or mixed s-zps-coordinate   ==! 
    … …  
    185197            ENDIF 
    186198            ! 
    187                DO jj = 1, jpj                         ! vertical interpolation of T & S 
     199            DO jj = 1, jpj                         ! vertical interpolation of T & S 
     200               DO ji = 1, jpi 
     201                  DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points 
     202                     zl = fsdept_n(ji,jj,jk) 
     203                     IF(     zl < gdept_1d(1  ) ) THEN         ! above the first level of data 
     204                        ztp(jk) = ztrcdta(ji,jj,1) 
     205                     ELSEIF( zl > gdept_1d(jpk) ) THEN         ! below the last level of data 
     206                        ztp(jk) =  ztrcdta(ji,jj,jpkm1) 
     207                     ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1 
     208                        DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1) 
     209                           IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN 
     210                              zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk)) 
     211                              ztp(jk) = ztrcdta(ji,jj,jkk) + ( ztrcdta(ji,jj,jkk+1) - & 
     212                                        ztrcdta(ji,jj,jkk) ) * zi  
     213                           ENDIF 
     214                        END DO 
     215                     ENDIF 
     216                  END DO 
     217                  DO jk = 1, jpkm1 
     218                    ztrcdta(ji,jj,jk) = ztp(jk) * tmask(ji,jj,jk)     ! mask required for mixed zps-s-coord 
     219                  END DO 
     220                  ztrcdta(ji,jj,jpk) = 0._wp 
     221                END DO 
     222            END DO 
     223            !  
     224         ELSE                                !==   z- or zps- coordinate   ==! 
     225            ! 
     226            IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level 
     227               DO jj = 1, jpj 
    188228                  DO ji = 1, jpi 
    189                      DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points 
    190                         zl = fsdept_n(ji,jj,jk) 
    191                         IF(     zl < gdept_1d(1  ) ) THEN         ! above the first level of data 
    192                            ztp(jk) =  sf_dta(1)%fnow(ji,jj,1) 
    193                         ELSEIF( zl > gdept_1d(jpk) ) THEN         ! below the last level of data 
    194                            ztp(jk) =  sf_dta(1)%fnow(ji,jj,jpkm1) 
    195                         ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1 
    196                            DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1) 
    197                               IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN 
    198                                  zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk)) 
    199                                  ztp(jk) = sf_dta(1)%fnow(ji,jj,jkk) + ( sf_dta(1)%fnow(ji,jj,jkk+1) - & 
    200                                            sf_dta(1)%fnow(ji,jj,jkk) ) * zi  
    201                               ENDIF 
    202                            END DO 
    203                         ENDIF 
    204                      END DO 
    205                      DO jk = 1, jpkm1 
    206                         sf_dta(1)%fnow(ji,jj,jk) = ztp(jk) * tmask(ji,jj,jk)     ! mask required for mixed zps-s-coord 
    207                      END DO 
    208                      sf_dta(1)%fnow(ji,jj,jpk) = 0._wp 
     229                     ik = mbkt(ji,jj)  
     230                     IF( ik > 1 ) THEN 
     231                        zl = ( gdept_1d(ik) - fsdept_n(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) ) 
     232                        ztrcdta(ji,jj,ik) = (1.-zl) * ztrcdta(ji,jj,ik) + zl * ztrcdta(ji,jj,ik-1) 
     233                     ENDIF 
     234                     ik = mikt(ji,jj) 
     235                     IF( ik > 1 ) THEN 
     236                        zl = ( fsdept_n(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) ) 
     237                        ztrcdta(ji,jj,ik) = (1.-zl) * ztrcdta(ji,jj,ik) + zl * ztrcdta(ji,jj,ik+1) 
     238                     ENDIF 
    209239                  END DO 
    210240               END DO 
    211             !  
    212          ELSE                                !==   z- or zps- coordinate   ==! 
    213             !                              
    214                sf_dta(1)%fnow(:,:,:) = sf_dta(1)%fnow(:,:,:) * tmask(:,:,:)    ! Mask 
    215                ! 
    216                IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level 
    217                   DO jj = 1, jpj 
    218                      DO ji = 1, jpi 
    219                         ik = mbkt(ji,jj)  
    220                         IF( ik > 1 ) THEN 
    221                            zl = ( gdept_1d(ik) - fsdept_n(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) ) 
    222                            sf_dta(1)%fnow(ji,jj,ik) = (1.-zl) * sf_dta(1)%fnow(ji,jj,ik) + zl * sf_dta(1)%fnow(ji,jj,ik-1) 
    223                         ENDIF 
    224                         ik = mikt(ji,jj) 
    225                         IF( ik > 1 ) THEN 
    226                            zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) ) 
    227                            sf_dta(1)%fnow(ji,jj,ik) = (1.-zl) * sf_dta(1)%fnow(ji,jj,ik) + zl * sf_dta(1)%fnow(ji,jj,ik+1) 
    228                         ENDIF 
    229                      END DO 
    230                   END DO 
    231                ENDIF 
    232             ! 
    233          ENDIF 
     241            ENDIF 
     242            ! 
     243         ENDIF 
     244         ! 
     245         ! Add multiplicative factor 
     246         ztrcdta(:,:,:) = ztrcdta(:,:,:) * ptrfac 
     247         ! 
     248         ! Data structure for trc_ini (and BFMv5.1 coupling) 
     249         IF( .NOT. PRESENT(ptrc) ) sf_dta(1)%fnow(:,:,:) = ztrcdta(:,:,:) 
     250         ! 
     251         ! Data structure for trc_dmp 
     252         IF( PRESENT(ptrc) )  ptrc(:,:,:) = ztrcdta(:,:,:) 
    234253         ! 
    235254         IF( lwp .AND. kt == nit000 ) THEN 
    … …  
    238257               WRITE(numout,*) 
    239258               WRITE(numout,*)'  level = 1' 
    240                CALL prihre( sf_dta(1)%fnow(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout ) 
     259               CALL prihre( ztrcdta(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout ) 
    241260               WRITE(numout,*)'  level = ', jpk/2 
    242                CALL prihre( sf_dta(1)%fnow(:,:,jpk/2), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout ) 
     261               CALL prihre( ztrcdta(:,:,jpk/2), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout ) 
    243262               WRITE(numout,*)'  level = ', jpkm1 
    244                CALL prihre( sf_dta(1)%fnow(:,:,jpkm1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout ) 
     263               CALL prihre( ztrcdta(:,:,jpkm1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout ) 
    245264               WRITE(numout,*) 
    246265         ENDIF 
     266         ! 
     267         CALL wrk_dealloc( jpi, jpj, jpk, ztrcdta ) 
     268         ! 
    247269      ENDIF 
    248270      ! 
    … …  
    255277   !!---------------------------------------------------------------------- 
    256278CONTAINS 
    257    SUBROUTINE trc_dta( kt, sf_dta, zrf_trfac )        ! Empty routine 
     279   SUBROUTINE trc_dta( kt, sf_dta, ptrfac, ptrc)        ! Empty routine 
    258280      WRITE(*,*) 'trc_dta: You should not have seen this print! error?', kt 
    259281   END SUBROUTINE trc_dta 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/TOP_SRC/trcini.F90

    r6308 r6746  
    123123               IF( ln_trc_ini(jn) ) THEN      ! update passive tracers arrays with input data read from file 
    124124                  jl = n_trc_index(jn)  
    125                   CALL trc_dta( nit000, sf_trcdta(jl) )   ! read tracer data at nit000 
    126                   trn(:,:,:,jn) = sf_trcdta(jl)%fnow(:,:,:) * tmask(:,:,:) * rf_trfac(jl) 
    127                   ! 
     125                  CALL trc_dta( nit000, sf_trcdta(jl), rf_trfac(jl) )   ! read tracer data at nit000 
     126                  trn(:,:,:,jn) = sf_trcdta(jl)%fnow(:,:,:)  
    128127                  IF( .NOT.ln_trcdmp .AND. .NOT.ln_trcdmp_clo ) THEN      !== deallocate data structure   ==! 
    129128                     !                                                    (data used only for initialisation) 
Note: See TracChangeset for help on using the changeset viewer.