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 6436 for branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO – NEMO

Ignore:
Timestamp:
2016-04-07T15:33:32+02:00 (8 years ago)
Author:
timgraham
Message:

Updated to r6424 of nemo_v3_6_STABLE

Location:
branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO
Files:
25 edited

Legend:

Unmodified
Added
Removed

  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r6333 r6436  
    253253   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fhld        !: heat flux from the lead used for bottom melting 
    254254 
    255    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw    !: snow-ocean mass exchange over 1 time step [kg/m2] 
    256    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_spr    !: snow precipitation on ice over 1 time step [kg/m2] 
    257    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sub    !: snow sublimation over 1 time step [kg/m2] 
    258  
    259    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_ice    !: ice-ocean mass exchange over 1 time step [kg/m2] 
    260    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sni    !: snow ice growth component of wfx_ice [kg/m2] 
    261    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_opw    !: lateral ice growth component of wfx_ice [kg/m2] 
    262    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_bog    !: bottom ice growth component of wfx_ice [kg/m2] 
    263    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_dyn    !: dynamical ice growth component of wfx_ice [kg/m2] 
    264    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_bom    !: bottom melt component of wfx_ice [kg/m2] 
    265    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sum    !: surface melt component of wfx_ice [kg/m2] 
    266    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_res    !: residual component of wfx_ice [kg/m2] 
    267  
    268    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_tot     !: ice concentration tendency (total) [s-1] 
     255   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw    !: snow-ocean mass exchange   [kg.m-2.s-1] 
     256   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_spr    !: snow precipitation on ice  [kg.m-2.s-1] 
     257   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sub    !: snow/ice sublimation       [kg.m-2.s-1] 
     258 
     259   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_ice    !: ice-ocean mass exchange                   [kg.m-2.s-1] 
     260   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sni    !: snow ice growth component of wfx_ice      [kg.m-2.s-1] 
     261   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_opw    !: lateral ice growth component of wfx_ice   [kg.m-2.s-1] 
     262   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_bog    !: bottom ice growth component of wfx_ice    [kg.m-2.s-1] 
     263   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_dyn    !: dynamical ice growth component of wfx_ice [kg.m-2.s-1] 
     264   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_bom    !: bottom melt component of wfx_ice          [kg.m-2.s-1] 
     265   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sum    !: surface melt component of wfx_ice         [kg.m-2.s-1] 
     266   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_res    !: residual component of wfx_ice             [kg.m-2.s-1] 
     267 
     268   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_tot     !: ice concentration tendency (total)          [s-1] 
    269269   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_thd     !: ice concentration tendency (thermodynamics) [s-1] 
    270    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_dyn     !: ice concentration tendency (dynamics) [s-1] 
     270   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_dyn     !: ice concentration tendency (dynamics)       [s-1] 
    271271 
    272272   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_bog     !: salt flux due to ice growth/melt                      [PSU/m2/s] 
    … …  
    279279   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_res     !: residual salt flux due to correction of ice thickness [PSU/m2/s] 
    280280 
    281    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_bog     !: total heat flux causing bottom ice growth  
    282    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_bom     !: total heat flux causing bottom ice melt  
    283    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_sum     !: total heat flux causing surface ice melt  
    284    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_opw     !: total heat flux causing open water ice formation 
    285    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_dif     !: total heat flux causing Temp change in the ice  
    286    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_snw     !: heat flux for snow melt  
    287    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err     !: heat flux error after heat diffusion  
    288    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_dif !: heat flux remaining due to change in non-solar flux 
    289    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_rem !: heat flux error after heat remapping  
    290    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_in      !: heat flux available for thermo transformations  
    291    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_out     !: heat flux remaining at the end of thermo transformations  
    292  
     281   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_sub     !: salt flux due to ice sublimation 
     282 
     283   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_bog     !: total heat flux causing bottom ice growth        [W.m-2] 
     284   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_bom     !: total heat flux causing bottom ice melt          [W.m-2] 
     285   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_sum     !: total heat flux causing surface ice melt         [W.m-2] 
     286   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_opw     !: total heat flux causing open water ice formation [W.m-2] 
     287   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_dif     !: total heat flux causing Temp change in the ice   [W.m-2] 
     288   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_snw     !: heat flux for snow melt                          [W.m-2] 
     289   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err     !: heat flux error after heat diffusion             [W.m-2] 
     290   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_dif !: heat flux remaining due to change in non-solar flux [W.m-2] 
     291   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_rem !: heat flux error after heat remapping             [W.m-2] 
     292   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_in      !: heat flux available for thermo transformations   [W.m-2] 
     293   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_out     !: heat flux remaining at the end of thermo transformations  [W.m-2] 
     294   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_err_sub !: mass flux error after sublimation [kg.m-2.s-1] 
     295    
    293296   ! heat flux associated with ice-atmosphere mass exchange 
    294    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_sub     !: heat flux for sublimation  
    295    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_spr     !: heat flux of the snow precipitation  
     297   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_sub     !: heat flux for sublimation  [W.m-2] 
     298   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_spr     !: heat flux of the snow precipitation  [W.m-2] 
    296299 
    297300   ! heat flux associated with ice-ocean mass exchange 
    298    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_thd     !: ice-ocean heat flux from thermo processes (limthd_dh)  
    299    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_dyn     !: ice-ocean heat flux from mecanical processes (limitd_me)  
    300    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_res     !: residual heat flux due to correction of ice thickness 
     301   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_thd     !: ice-ocean heat flux from thermo processes (limthd_dh)  [W.m-2] 
     302   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_dyn     !: ice-ocean heat flux from mecanical processes (limitd_me)  [W.m-2] 
     303   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_res     !: residual heat flux due to correction of ice thickness [W.m-2] 
    301304 
    302305   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ftr_ice   !: transmitted solar radiation under ice 
    … …  
    441444         &      fhtur  (jpi,jpj) , ftr_ice(jpi,jpj,jpl), qlead  (jpi,jpj) ,                     & 
    442445         &      rn_amax_2d(jpi,jpj),                                                            & 
    443          &      sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) ,                        & 
     446         &      sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) , sfx_sub(jpi,jpj) ,                       & 
    444447         &      sfx_bog(jpi,jpj) , sfx_bom(jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) ,    & 
    445448         &      hfx_res(jpi,jpj) , hfx_snw(jpi,jpj) , hfx_sub(jpi,jpj) , hfx_err(jpi,jpj) ,     &  
    446          &      hfx_err_dif(jpi,jpj) , hfx_err_rem(jpi,jpj) ,                                   & 
     449         &      hfx_err_dif(jpi,jpj) , hfx_err_rem(jpi,jpj) , wfx_err_sub(jpi,jpj) ,       & 
    447450         &      hfx_in (jpi,jpj) , hfx_out(jpi,jpj) , fhld(jpi,jpj) ,                           & 
    448451         &      hfx_sum(jpi,jpj) , hfx_bom(jpi,jpj) , hfx_bog(jpi,jpj) , hfx_dif(jpi,jpj) , hfx_opw(jpi,jpj) ,    & 

  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90

    r6333 r6436  
    2424   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    2525   USE sbc_oce , ONLY : sfx  ! Surface boundary condition: ocean fields 
    26  
     26   USE sbc_ice , ONLY : qevap_ice 
     27    
    2728   IMPLICIT NONE 
    2829   PRIVATE 
    … …  
    184185         ! salt flux 
    185186         zfs_b  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
    186             &                  sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:)                                  & 
     187            &                  sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:)                   & 
    187188            &                ) *  e12t(:,:) * tmask(:,:,1) * zconv ) 
    188189 
    … …  
    209210         ! salt flux 
    210211         zfs  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
    211             &                sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:)                                  &  
     212            &                sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:)                   &  
    212213            &              ) * e12t(:,:) * tmask(:,:,1) * zconv ) - zfs_b 
    213214 
    … …  
    287288#if ! defined key_bdy 
    288289      ! heat flux 
    289       zhfx  = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es - hfx_sub ) * e12t * tmask(:,:,1) * zconv )  
     290      zhfx  = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es - SUM( qevap_ice * a_i_b, dim=3 ) )  & 
     291         &              * e12t * tmask(:,:,1) * zconv )  
    290292      ! salt flux 
    291293      zsfx  = glob_sum( ( sfx + diag_smvi ) * e12t * tmask(:,:,1) * zconv ) * rday 

  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90

    r5781 r6436  
    5656      real(wp)   ::   zbg_ivo, zbg_svo, zbg_are, zbg_sal ,zbg_tem ,zbg_ihc ,zbg_shc 
    5757      real(wp)   ::   zbg_sfx, zbg_sfx_bri, zbg_sfx_bog, zbg_sfx_bom, zbg_sfx_sum, zbg_sfx_sni,   & 
    58       &               zbg_sfx_opw, zbg_sfx_res, zbg_sfx_dyn  
     58      &               zbg_sfx_opw, zbg_sfx_res, zbg_sfx_dyn, zbg_sfx_sub  
    5959      real(wp)   ::   zbg_vfx, zbg_vfx_bog, zbg_vfx_opw, zbg_vfx_sni, zbg_vfx_dyn 
    6060      real(wp)   ::   zbg_vfx_bom, zbg_vfx_sum, zbg_vfx_res, zbg_vfx_spr, zbg_vfx_snw, zbg_vfx_sub   
    … …  
    111111      zbg_sfx_bom = ztmp * glob_sum( sfx_bom(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    112112      zbg_sfx_sum = ztmp * glob_sum( sfx_sum(:,:) * e12t(:,:) * tmask(:,:,1) ) 
     113      zbg_sfx_sub = ztmp * glob_sum( sfx_sub(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    113114 
    114115      ! Heat budget 
    … …  
    189190      CALL iom_put( 'ibgsfxbom' , zbg_sfx_bom                              )   ! salt flux bottom melt       - 
    190191      CALL iom_put( 'ibgsfxsum' , zbg_sfx_sum                              )   ! salt flux surface melt      - 
     192      CALL iom_put( 'ibgsfxsub' , zbg_sfx_sub                              )   ! salt flux sublimation      - 
    191193 
    192194      CALL iom_put( 'ibghfxdhc' , zbg_hfx_dhc                              )   ! Heat content variation in snow and ice [W] 

  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90

    r5781 r6436  
    4545   REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   asum     ! sum of total ice and open water area 
    4646   REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   aksum    ! ratio of area removed to area ridged 
    47    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   athorn   ! participation function; fraction of ridging/ 
    48    !                                                     ! closing associated w/ category n 
     47   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   athorn   ! participation function; fraction of ridging/closing associated w/ category n 
    4948   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hrmin    ! minimum ridge thickness 
    5049   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hrmax    ! maximum ridge thickness 
    5150   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hraft    ! thickness of rafted ice 
    52    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   krdg     ! mean ridge thickness/thickness of ridging ice  
     51   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   krdg     ! thickness of ridging ice / mean ridge thickness 
    5352   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   aridge   ! participating ice ridging 
    5453   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   araft    ! participating ice rafting 
    5554 
    5655   REAL(wp), PARAMETER ::   krdgmin = 1.1_wp    ! min ridge thickness multiplier 
    57    REAL(wp), PARAMETER ::   kraft   = 2.0_wp    ! rafting multipliyer 
    58    REAL(wp), PARAMETER ::   kamax   = 1.0_wp    ! max of ice area authorized (clem: scheme is not stable if kamax <= 0.99) 
     56   REAL(wp), PARAMETER ::   kraft   = 0.5_wp    ! rafting multipliyer 
    5957 
    6058   REAL(wp) ::   Cp                             !  
    6159   ! 
    62    !----------------------------------------------------------------------- 
    63    ! Ridging diagnostic arrays for history files 
    64    !----------------------------------------------------------------------- 
    65    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dardg1dt   ! rate of fractional area loss by ridging ice (1/s) 
    66    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dardg2dt   ! rate of fractional area gain by new ridges (1/s) 
    67    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dvirdgdt   ! rate of ice volume ridged (m/s) 
    68    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   opening    ! rate of opening due to divergence/shear (1/s) 
    6960   ! 
    7061   !!---------------------------------------------------------------------- 
    … …  
    8374         &      asum (jpi,jpj)     , athorn(jpi,jpj,0:jpl)                    ,     & 
    8475         &      aksum(jpi,jpj)                                                ,     & 
    85          ! 
    8676         &      hrmin(jpi,jpj,jpl) , hraft(jpi,jpj,jpl) , aridge(jpi,jpj,jpl) ,     & 
    87          &      hrmax(jpi,jpj,jpl) , krdg (jpi,jpj,jpl) , araft (jpi,jpj,jpl) ,     & 
    88          ! 
    89          !* Ridging diagnostic arrays for history files 
    90          &      dardg1dt(jpi,jpj)  , dardg2dt(jpi,jpj)                        ,     &  
    91          &      dvirdgdt(jpi,jpj)  , opening(jpi,jpj)                         , STAT=lim_itd_me_alloc ) 
     77         &      hrmax(jpi,jpj,jpl) , krdg (jpi,jpj,jpl) , araft (jpi,jpj,jpl) , STAT=lim_itd_me_alloc ) 
    9278         ! 
    9379      IF( lim_itd_me_alloc /= 0 )   CALL ctl_warn( 'lim_itd_me_alloc: failed to allocate arrays' ) 
    … …  
    132118      REAL(wp), POINTER, DIMENSION(:,:)   ::   opning          ! rate of opening due to divergence/shear 
    133119      REAL(wp), POINTER, DIMENSION(:,:)   ::   closing_gross   ! rate at which area removed, not counting area of new ridges 
    134       REAL(wp), POINTER, DIMENSION(:,:)   ::   msnow_mlt       ! mass of snow added to ocean (kg m-2) 
    135       REAL(wp), POINTER, DIMENSION(:,:)   ::   esnow_mlt       ! energy needed to melt snow in ocean (J m-2) 
    136       REAL(wp), POINTER, DIMENSION(:,:)   ::   vt_i_init, vt_i_final  !  ice volume summed over categories 
    137120      ! 
    138121      INTEGER, PARAMETER ::   nitermax = 20     
    … …  
    142125      IF( nn_timing == 1 )  CALL timing_start('limitd_me') 
    143126 
    144       CALL wrk_alloc( jpi,jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 
     127      CALL wrk_alloc( jpi,jpj, closing_net, divu_adv, opning, closing_gross ) 
    145128 
    146129      IF(ln_ctl) THEN 
    … …  
    154137      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    155138 
    156       CALL lim_var_zapsmall 
    157       CALL lim_var_glo2eqv            ! equivalent variables, requested for rafting 
    158  
    159139      !-----------------------------------------------------------------------------! 
    160140      ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons 
    … …  
    164144      CALL lim_itd_me_ridgeprep                                    ! prepare ridging 
    165145      ! 
    166       IF( con_i)   CALL lim_column_sum( jpl, v_i, vt_i_init )      ! conservation check 
    167146 
    168147      DO jj = 1, jpj                                     ! Initialize arrays. 
    169148         DO ji = 1, jpi 
    170             msnow_mlt(ji,jj) = 0._wp 
    171             esnow_mlt(ji,jj) = 0._wp 
    172             dardg1dt (ji,jj) = 0._wp 
    173             dardg2dt (ji,jj) = 0._wp 
    174             dvirdgdt (ji,jj) = 0._wp 
    175             opening  (ji,jj) = 0._wp 
    176149 
    177150            !-----------------------------------------------------------------------------! 
    … …  
    204177            ! If divu_adv < 0, make sure the closing rate is large enough 
    205178            ! to give asum = 1.0 after ridging. 
    206  
    207             divu_adv(ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice  ! asum found in ridgeprep 
     179             
     180            divu_adv(ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice  ! asum found in ridgeprep 
    208181 
    209182            IF( divu_adv(ji,jj) < 0._wp )   closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) ) 
    … …  
    224197      DO WHILE ( iterate_ridging > 0 .AND. niter < nitermax ) 
    225198 
     199         ! 3.2 closing_gross 
     200         !-----------------------------------------------------------------------------! 
     201         ! Based on the ITD of ridging and ridged ice, convert the net 
     202         !  closing rate to a gross closing rate.   
     203         ! NOTE: 0 < aksum <= 1 
     204         closing_gross(:,:) = closing_net(:,:) / aksum(:,:) 
     205 
     206         ! correction to closing rate and opening if closing rate is excessive 
     207         !--------------------------------------------------------------------- 
     208         ! Reduce the closing rate if more than 100% of the open water  
     209         ! would be removed.  Reduce the opening rate proportionately. 
    226210         DO jj = 1, jpj 
    227211            DO ji = 1, jpi 
    228  
    229                ! 3.2 closing_gross 
    230                !-----------------------------------------------------------------------------! 
    231                ! Based on the ITD of ridging and ridged ice, convert the net 
    232                !  closing rate to a gross closing rate.   
    233                ! NOTE: 0 < aksum <= 1 
    234                closing_gross(ji,jj) = closing_net(ji,jj) / aksum(ji,jj) 
    235  
    236                ! correction to closing rate and opening if closing rate is excessive 
    237                !--------------------------------------------------------------------- 
    238                ! Reduce the closing rate if more than 100% of the open water  
    239                ! would be removed.  Reduce the opening rate proportionately. 
    240                za   = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice 
    241                IF( za > epsi20 ) THEN 
    242                   zfac = MIN( 1._wp, ato_i(ji,jj) / za ) 
    243                   closing_gross(ji,jj) = closing_gross(ji,jj) * zfac 
    244                   opning       (ji,jj) = opning       (ji,jj) * zfac 
     212               za   = ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice 
     213               IF( za < 0._wp .AND. za > - ato_i(ji,jj) ) THEN  ! would lead to negative ato_i 
     214                  zfac = - ato_i(ji,jj) / za 
     215                  opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) - ato_i(ji,jj) * r1_rdtice  
     216               ELSEIF( za > 0._wp .AND. za > ( asum(ji,jj) - ato_i(ji,jj) ) ) THEN  ! would lead to ato_i > asum 
     217                  zfac = ( asum(ji,jj) - ato_i(ji,jj) ) / za 
     218                  opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) + ( asum(ji,jj) - ato_i(ji,jj) ) * r1_rdtice  
    245219               ENDIF 
    246  
    247220            END DO 
    248221         END DO 
    … …  
    256229               DO ji = 1, jpi 
    257230                  za = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice 
    258                   IF( za  >  epsi20 ) THEN 
    259                      zfac = MIN( 1._wp, a_i(ji,jj,jl) / za ) 
     231                  IF( za  >  a_i(ji,jj,jl) ) THEN 
     232                     zfac = a_i(ji,jj,jl) / za 
    260233                     closing_gross(ji,jj) = closing_gross(ji,jj) * zfac 
    261                      opning       (ji,jj) = opning       (ji,jj) * zfac 
    262234                  ENDIF 
    263235               END DO 
    … …  
    268240         !-----------------------------------------------------------------------------! 
    269241 
    270          CALL lim_itd_me_ridgeshift( opning, closing_gross, msnow_mlt, esnow_mlt ) 
    271  
     242         CALL lim_itd_me_ridgeshift( opning, closing_gross ) 
     243 
     244          
    272245         ! 3.4 Compute total area of ice plus open water after ridging. 
    273246         !-----------------------------------------------------------------------------! 
    274247         ! This is in general not equal to one because of divergence during transport 
    275          asum(:,:) = ato_i(:,:) 
    276          DO jl = 1, jpl 
    277             asum(:,:) = asum(:,:) + a_i(:,:,jl) 
    278          END DO 
     248         asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 ) 
    279249 
    280250         ! 3.5 Do we keep on iterating ??? 
    … …  
    284254 
    285255         iterate_ridging = 0 
    286  
    287256         DO jj = 1, jpj 
    288257            DO ji = 1, jpi 
    289                IF (ABS(asum(ji,jj) - kamax ) < epsi10) THEN 
     258               IF( ABS( asum(ji,jj) - 1._wp ) < epsi10 ) THEN 
    290259                  closing_net(ji,jj) = 0._wp 
    291260                  opning     (ji,jj) = 0._wp 
    292261               ELSE 
    293262                  iterate_ridging    = 1 
    294                   divu_adv   (ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice 
     263                  divu_adv   (ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice 
    295264                  closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) ) 
    296265                  opning     (ji,jj) = MAX( 0._wp,  divu_adv(ji,jj) ) 
    … …  
    309278 
    310279         IF( iterate_ridging == 1 ) THEN 
     280            CALL lim_itd_me_ridgeprep 
    311281            IF( niter  >  nitermax ) THEN 
    312282               WRITE(numout,*) ' ALERTE : non-converging ridging scheme ' 
    313283               WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging 
    314284            ENDIF 
    315             CALL lim_itd_me_ridgeprep 
    316285         ENDIF 
    317286 
    318287      END DO !! on the do while over iter 
    319  
    320       !-----------------------------------------------------------------------------! 
    321       ! 4) Ridging diagnostics 
    322       !-----------------------------------------------------------------------------! 
    323       ! Convert ridging rate diagnostics to correct units. 
    324       ! Update fresh water and heat fluxes due to snow melt. 
    325       DO jj = 1, jpj 
    326          DO ji = 1, jpi 
    327  
    328             dardg1dt(ji,jj) = dardg1dt(ji,jj) * r1_rdtice 
    329             dardg2dt(ji,jj) = dardg2dt(ji,jj) * r1_rdtice 
    330             dvirdgdt(ji,jj) = dvirdgdt(ji,jj) * r1_rdtice 
    331             opening (ji,jj) = opening (ji,jj) * r1_rdtice 
    332  
    333             !-----------------------------------------------------------------------------! 
    334             ! 5) Heat, salt and freshwater fluxes 
    335             !-----------------------------------------------------------------------------! 
    336             wfx_snw(ji,jj) = wfx_snw(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice     ! fresh water source for ocean 
    337             hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice     ! heat sink for ocean (<0, W.m-2) 
    338  
    339          END DO 
    340       END DO 
    341  
    342       ! Check if there is a ridging error 
    343       IF( lwp ) THEN 
    344          DO jj = 1, jpj 
    345             DO ji = 1, jpi 
    346                IF( ABS( asum(ji,jj) - kamax)  >  epsi10 ) THEN   ! there is a bug 
    347                   WRITE(numout,*) ' ' 
    348                   WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj) 
    349                   WRITE(numout,*) ' limitd_me ' 
    350                   WRITE(numout,*) ' POINT : ', ji, jj 
    351                   WRITE(numout,*) ' jpl, a_i, athorn ' 
    352                   WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0) 
    353                   DO jl = 1, jpl 
    354                      WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl) 
    355                   END DO 
    356                ENDIF 
    357             END DO 
    358          END DO 
    359       END IF 
    360  
    361       ! Conservation check 
    362       IF ( con_i ) THEN 
    363          CALL lim_column_sum (jpl,   v_i, vt_i_final) 
    364          fieldid = ' v_i : limitd_me ' 
    365          CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)  
    366       ENDIF 
    367288 
    368289      CALL lim_var_agg( 1 )  
    … …  
    410331      ENDIF  ! ln_limdyn=.true. 
    411332      ! 
    412       CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 
     333      CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross ) 
    413334      ! 
    414335      IF( nn_timing == 1 )  CALL timing_stop('limitd_me') 
    415336   END SUBROUTINE lim_itd_me 
    416337 
     338   SUBROUTINE lim_itd_me_ridgeprep 
     339      !!---------------------------------------------------------------------! 
     340      !!                ***  ROUTINE lim_itd_me_ridgeprep *** 
     341      !! 
     342      !! ** Purpose :   preparation for ridging and strength calculations 
     343      !! 
     344      !! ** Method  :   Compute the thickness distribution of the ice and open water  
     345      !!              participating in ridging and of the resulting ridges. 
     346      !!---------------------------------------------------------------------! 
     347      INTEGER ::   ji,jj, jl    ! dummy loop indices 
     348      REAL(wp) ::   Gstari, astari, hrmean, zdummy   ! local scalar 
     349      REAL(wp), POINTER, DIMENSION(:,:,:) ::   Gsum      ! Gsum(n) = sum of areas in categories 0 to n 
     350      !------------------------------------------------------------------------------! 
     351 
     352      CALL wrk_alloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 ) 
     353 
     354      Gstari     = 1.0/rn_gstar     
     355      astari     = 1.0/rn_astar     
     356      aksum(:,:)    = 0.0 
     357      athorn(:,:,:) = 0.0 
     358      aridge(:,:,:) = 0.0 
     359      araft (:,:,:) = 0.0 
     360 
     361      ! Zero out categories with very small areas 
     362      CALL lim_var_zapsmall 
     363 
     364      ! Ice thickness needed for rafting 
     365      DO jl = 1, jpl 
     366         DO jj = 1, jpj 
     367            DO ji = 1, jpi 
     368               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) 
     369               ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 
     370            END DO 
     371         END DO 
     372      END DO 
     373 
     374      !------------------------------------------------------------------------------! 
     375      ! 1) Participation function  
     376      !------------------------------------------------------------------------------! 
     377 
     378      ! Compute total area of ice plus open water. 
     379      ! This is in general not equal to one because of divergence during transport 
     380      asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 ) 
     381 
     382      ! Compute cumulative thickness distribution function 
     383      ! Compute the cumulative thickness distribution function Gsum, 
     384      ! where Gsum(n) is the fractional area in categories 0 to n. 
     385      ! initial value (in h = 0) equals open water area 
     386      Gsum(:,:,-1) = 0._wp 
     387      Gsum(:,:,0 ) = ato_i(:,:) 
     388      ! for each value of h, you have to add ice concentration then 
     389      DO jl = 1, jpl 
     390         Gsum(:,:,jl) = Gsum(:,:,jl-1) + a_i(:,:,jl) 
     391      END DO 
     392 
     393      ! Normalize the cumulative distribution to 1 
     394      DO jl = 0, jpl 
     395         Gsum(:,:,jl) = Gsum(:,:,jl) / asum(:,:) 
     396      END DO 
     397 
     398      ! 1.3 Compute participation function a(h) = b(h).g(h) (athorn) 
     399      !-------------------------------------------------------------------------------------------------- 
     400      ! Compute the participation function athorn; this is analogous to 
     401      ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975). 
     402      ! area lost from category n due to ridging/closing 
     403      ! athorn(n)   = total area lost due to ridging/closing 
     404      ! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).  
     405      ! 
     406      ! The expressions for athorn are found by integrating b(h)g(h) between 
     407      ! the category boundaries. 
     408      ! athorn is always >= 0 and SUM(athorn(0:jpl))=1 
     409      !----------------------------------------------------------------- 
     410 
     411      IF( nn_partfun == 0 ) THEN       !--- Linear formulation (Thorndike et al., 1975) 
     412         DO jl = 0, jpl     
     413            DO jj = 1, jpj  
     414               DO ji = 1, jpi 
     415                  IF    ( Gsum(ji,jj,jl)   < rn_gstar ) THEN 
     416                     athorn(ji,jj,jl) = Gstari * ( Gsum(ji,jj,jl) - Gsum(ji,jj,jl-1) ) * & 
     417                        &                        ( 2._wp - ( Gsum(ji,jj,jl-1) + Gsum(ji,jj,jl) ) * Gstari ) 
     418                  ELSEIF( Gsum(ji,jj,jl-1) < rn_gstar ) THEN 
     419                     athorn(ji,jj,jl) = Gstari * ( rn_gstar       - Gsum(ji,jj,jl-1) ) *  & 
     420                        &                        ( 2._wp - ( Gsum(ji,jj,jl-1) + rn_gstar       ) * Gstari ) 
     421                  ELSE 
     422                     athorn(ji,jj,jl) = 0._wp 
     423                  ENDIF 
     424               END DO 
     425            END DO 
     426         END DO 
     427 
     428      ELSE                             !--- Exponential, more stable formulation (Lipscomb et al, 2007) 
     429         !                         
     430         zdummy = 1._wp / ( 1._wp - EXP(-astari) )        ! precompute exponential terms using Gsum as a work array 
     431         DO jl = -1, jpl 
     432            Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy 
     433         END DO 
     434         DO jl = 0, jpl 
     435             athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl) 
     436         END DO 
     437         ! 
     438      ENDIF 
     439 
     440      IF( ln_rafting ) THEN      ! Ridging and rafting ice participation functions 
     441         ! 
     442         DO jl = 1, jpl 
     443            DO jj = 1, jpj  
     444               DO ji = 1, jpi 
     445                  zdummy           = TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) 
     446                  aridge(ji,jj,jl) = ( 1._wp + zdummy ) * 0.5_wp * athorn(ji,jj,jl) 
     447                  araft (ji,jj,jl) = ( 1._wp - zdummy ) * 0.5_wp * athorn(ji,jj,jl) 
     448               END DO 
     449            END DO 
     450         END DO 
     451 
     452      ELSE 
     453         ! 
     454         DO jl = 1, jpl 
     455            aridge(:,:,jl) = athorn(:,:,jl) 
     456         END DO 
     457         ! 
     458      ENDIF 
     459 
     460      !----------------------------------------------------------------- 
     461      ! 2) Transfer function 
     462      !----------------------------------------------------------------- 
     463      ! Compute max and min ridged ice thickness for each ridging category. 
     464      ! Assume ridged ice is uniformly distributed between hrmin and hrmax. 
     465      !  
     466      ! This parameterization is a modified version of Hibler (1980). 
     467      ! The mean ridging thickness, hrmean, is proportional to hi^(0.5) 
     468      !  and for very thick ridging ice must be >= krdgmin*hi 
     469      ! 
     470      ! The minimum ridging thickness, hrmin, is equal to 2*hi  
     471      !  (i.e., rafting) and for very thick ridging ice is 
     472      !  constrained by hrmin <= (hrmean + hi)/2. 
     473      !  
     474      ! The maximum ridging thickness, hrmax, is determined by 
     475      !  hrmean and hrmin. 
     476      ! 
     477      ! These modifications have the effect of reducing the ice strength 
     478      ! (relative to the Hibler formulation) when very thick ice is 
     479      ! ridging. 
     480      ! 
     481      ! aksum = net area removed/ total area removed 
     482      ! where total area removed = area of ice that ridges 
     483      !         net area removed = total area removed - area of new ridges 
     484      !----------------------------------------------------------------- 
     485 
     486      aksum(:,:) = athorn(:,:,0) 
     487      ! Transfer function 
     488      DO jl = 1, jpl !all categories have a specific transfer function 
     489         DO jj = 1, jpj 
     490            DO ji = 1, jpi 
     491                
     492               IF( athorn(ji,jj,jl) > 0._wp ) THEN 
     493                  hrmean          = MAX( SQRT( rn_hstar * ht_i(ji,jj,jl) ), ht_i(ji,jj,jl) * krdgmin ) 
     494                  hrmin(ji,jj,jl) = MIN( 2._wp * ht_i(ji,jj,jl), 0.5_wp * ( hrmean + ht_i(ji,jj,jl) ) ) 
     495                  hrmax(ji,jj,jl) = 2._wp * hrmean - hrmin(ji,jj,jl) 
     496                  hraft(ji,jj,jl) = ht_i(ji,jj,jl) / kraft 
     497                  krdg(ji,jj,jl)  = ht_i(ji,jj,jl) / MAX( hrmean, epsi20 ) 
     498 
     499                  ! Normalization factor : aksum, ensures mass conservation 
     500                  aksum(ji,jj) = aksum(ji,jj) + aridge(ji,jj,jl) * ( 1._wp - krdg(ji,jj,jl) )    & 
     501                     &                        + araft (ji,jj,jl) * ( 1._wp - kraft          ) 
     502 
     503               ELSE 
     504                  hrmin(ji,jj,jl)  = 0._wp  
     505                  hrmax(ji,jj,jl)  = 0._wp  
     506                  hraft(ji,jj,jl)  = 0._wp  
     507                  krdg (ji,jj,jl)  = 1._wp 
     508               ENDIF 
     509 
     510            END DO 
     511         END DO 
     512      END DO 
     513      ! 
     514      CALL wrk_dealloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 ) 
     515      ! 
     516   END SUBROUTINE lim_itd_me_ridgeprep 
     517 
     518 
     519   SUBROUTINE lim_itd_me_ridgeshift( opning, closing_gross ) 
     520      !!---------------------------------------------------------------------- 
     521      !!                ***  ROUTINE lim_itd_me_icestrength *** 
     522      !! 
     523      !! ** Purpose :   shift ridging ice among thickness categories of ice thickness 
     524      !! 
     525      !! ** Method  :   Remove area, volume, and energy from each ridging category 
     526      !!              and add to thicker ice categories. 
     527      !!---------------------------------------------------------------------- 
     528      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   opning         ! rate of opening due to divergence/shear 
     529      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   closing_gross  ! rate at which area removed, excluding area of new ridges 
     530      ! 
     531      CHARACTER (len=80) ::   fieldid   ! field identifier 
     532      ! 
     533      INTEGER ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices 
     534      INTEGER ::   ij                ! horizontal index, combines i and j loops 
     535      INTEGER ::   icells            ! number of cells with a_i > puny 
     536      REAL(wp) ::   hL, hR, farea    ! left and right limits of integration 
     537 
     538      INTEGER , POINTER, DIMENSION(:) ::   indxi, indxj   ! compressed indices 
     539      REAL(wp), POINTER, DIMENSION(:) ::   zswitch, fvol   ! new ridge volume going to n2 
     540 
     541      REAL(wp), POINTER, DIMENSION(:) ::   afrac            ! fraction of category area ridged  
     542      REAL(wp), POINTER, DIMENSION(:) ::   ardg1 , ardg2    ! area of ice ridged & new ridges 
     543      REAL(wp), POINTER, DIMENSION(:) ::   vsrdg , esrdg    ! snow volume & energy of ridging ice 
     544      REAL(wp), POINTER, DIMENSION(:) ::   dhr   , dhr2     ! hrmax - hrmin  &  hrmax^2 - hrmin^2 
     545 
     546      REAL(wp), POINTER, DIMENSION(:) ::   vrdg1   ! volume of ice ridged 
     547      REAL(wp), POINTER, DIMENSION(:) ::   vrdg2   ! volume of new ridges 
     548      REAL(wp), POINTER, DIMENSION(:) ::   vsw     ! volume of seawater trapped into ridges 
     549      REAL(wp), POINTER, DIMENSION(:) ::   srdg1   ! sal*volume of ice ridged 
     550      REAL(wp), POINTER, DIMENSION(:) ::   srdg2   ! sal*volume of new ridges 
     551      REAL(wp), POINTER, DIMENSION(:) ::   smsw    ! sal*volume of water trapped into ridges 
     552      REAL(wp), POINTER, DIMENSION(:) ::   oirdg1, oirdg2   ! ice age of ice ridged 
     553 
     554      REAL(wp), POINTER, DIMENSION(:) ::   afrft            ! fraction of category area rafted 
     555      REAL(wp), POINTER, DIMENSION(:) ::   arft1 , arft2    ! area of ice rafted and new rafted zone 
     556      REAL(wp), POINTER, DIMENSION(:) ::   virft , vsrft    ! ice & snow volume of rafting ice 
     557      REAL(wp), POINTER, DIMENSION(:) ::   esrft , smrft    ! snow energy & salinity of rafting ice 
     558      REAL(wp), POINTER, DIMENSION(:) ::   oirft1, oirft2   ! ice age of ice rafted 
     559 
     560      REAL(wp), POINTER, DIMENSION(:,:) ::   eirft      ! ice energy of rafting ice 
     561      REAL(wp), POINTER, DIMENSION(:,:) ::   erdg1      ! enth*volume of ice ridged 
     562      REAL(wp), POINTER, DIMENSION(:,:) ::   erdg2      ! enth*volume of new ridges 
     563      REAL(wp), POINTER, DIMENSION(:,:) ::   ersw       ! enth of water trapped into ridges 
     564      !!---------------------------------------------------------------------- 
     565 
     566      CALL wrk_alloc( jpij,        indxi, indxj ) 
     567      CALL wrk_alloc( jpij,        zswitch, fvol ) 
     568      CALL wrk_alloc( jpij,        afrac, ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 ) 
     569      CALL wrk_alloc( jpij,        vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw, oirdg1, oirdg2 ) 
     570      CALL wrk_alloc( jpij,        afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 ) 
     571      CALL wrk_alloc( jpij,nlay_i, eirft, erdg1, erdg2, ersw ) 
     572 
     573      !------------------------------------------------------------------------------- 
     574      ! 1) Compute change in open water area due to closing and opening. 
     575      !------------------------------------------------------------------------------- 
     576      DO jj = 1, jpj 
     577         DO ji = 1, jpi 
     578            ato_i(ji,jj) = MAX( 0._wp, ato_i(ji,jj) +  & 
     579               &                     ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice ) 
     580         END DO 
     581      END DO 
     582 
     583      !----------------------------------------------------------------- 
     584      ! 3) Pump everything from ice which is being ridged / rafted 
     585      !----------------------------------------------------------------- 
     586      ! Compute the area, volume, and energy of ice ridging in each 
     587      ! category, along with the area of the resulting ridge. 
     588 
     589      DO jl1 = 1, jpl !jl1 describes the ridging category 
     590 
     591         !------------------------------------------------ 
     592         ! 3.1) Identify grid cells with nonzero ridging 
     593         !------------------------------------------------ 
     594         icells = 0 
     595         DO jj = 1, jpj 
     596            DO ji = 1, jpi 
     597               IF( athorn(ji,jj,jl1) > 0._wp .AND. closing_gross(ji,jj) > 0._wp ) THEN 
     598                  icells = icells + 1 
     599                  indxi(icells) = ji 
     600                  indxj(icells) = jj 
     601               ENDIF 
     602            END DO 
     603         END DO 
     604 
     605         DO ij = 1, icells 
     606            ji = indxi(ij) ; jj = indxj(ij) 
     607 
     608            !-------------------------------------------------------------------- 
     609            ! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2) 
     610            !-------------------------------------------------------------------- 
     611            ardg1(ij) = aridge(ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice 
     612            arft1(ij) = araft (ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice 
     613 
     614            !--------------------------------------------------------------- 
     615            ! 3.3) Compute ridging /rafting fractions, make sure afrac <=1  
     616            !--------------------------------------------------------------- 
     617            afrac(ij) = ardg1(ij) / a_i(ji,jj,jl1) !ridging 
     618            afrft(ij) = arft1(ij) / a_i(ji,jj,jl1) !rafting 
     619            ardg2(ij) = ardg1(ij) * krdg(ji,jj,jl1) 
     620            arft2(ij) = arft1(ij) * kraft 
     621 
     622            !-------------------------------------------------------------------------- 
     623            ! 3.4) Subtract area, volume, and energy from ridging  
     624            !     / rafting category n1. 
     625            !-------------------------------------------------------------------------- 
     626            vrdg1(ij) = v_i(ji,jj,jl1) * afrac(ij) 
     627            vrdg2(ij) = vrdg1(ij) * ( 1. + rn_por_rdg ) 
     628            vsw  (ij) = vrdg1(ij) * rn_por_rdg 
     629 
     630            vsrdg (ij) = v_s  (ji,jj,  jl1) * afrac(ij) 
     631            esrdg (ij) = e_s  (ji,jj,1,jl1) * afrac(ij) 
     632            srdg1 (ij) = smv_i(ji,jj,  jl1) * afrac(ij) 
     633            oirdg1(ij) = oa_i (ji,jj,  jl1) * afrac(ij) 
     634            oirdg2(ij) = oa_i (ji,jj,  jl1) * afrac(ij) * krdg(ji,jj,jl1)  
     635 
     636            ! rafting volumes, heat contents ... 
     637            virft (ij) = v_i  (ji,jj,  jl1) * afrft(ij) 
     638            vsrft (ij) = v_s  (ji,jj,  jl1) * afrft(ij) 
     639            esrft (ij) = e_s  (ji,jj,1,jl1) * afrft(ij) 
     640            smrft (ij) = smv_i(ji,jj,  jl1) * afrft(ij)  
     641            oirft1(ij) = oa_i (ji,jj,  jl1) * afrft(ij)  
     642            oirft2(ij) = oa_i (ji,jj,  jl1) * afrft(ij) * kraft  
     643 
     644            !----------------------------------------------------------------- 
     645            ! 3.5) Compute properties of new ridges 
     646            !----------------------------------------------------------------- 
     647            smsw(ij)  = vsw(ij) * sss_m(ji,jj)                   ! salt content of seawater frozen in voids 
     648            srdg2(ij) = srdg1(ij) + smsw(ij)                     ! salt content of new ridge 
     649             
     650            sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ij) * rhoic * r1_rdtice 
     651            wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ij) * rhoic * r1_rdtice   ! increase in ice volume due to seawater frozen in voids 
     652             
     653            !------------------------------------------             
     654            ! 3.7 Put the snow somewhere in the ocean 
     655            !------------------------------------------             
     656            !  Place part of the snow lost by ridging into the ocean.  
     657            !  Note that esrdg > 0; the ocean must cool to melt snow. 
     658            !  If the ocean temp = Tf already, new ice must grow. 
     659            !  During the next time step, thermo_rates will determine whether 
     660            !  the ocean cools or new ice grows. 
     661            wfx_snw(ji,jj) = wfx_snw(ji,jj) + ( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnowrdg )   &  
     662               &                              + rhosn * vsrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice  ! fresh water source for ocean 
     663 
     664            hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ( - esrdg(ij) * ( 1._wp - rn_fsnowrdg )         &  
     665               &                                - esrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice        ! heat sink for ocean (<0, W.m-2) 
     666 
     667            !----------------------------------------------------------------- 
     668            ! 3.8 Compute quantities used to apportion ice among categories 
     669            ! in the n2 loop below 
     670            !----------------------------------------------------------------- 
     671            dhr (ij) = 1._wp / ( hrmax(ji,jj,jl1)                    - hrmin(ji,jj,jl1)                    ) 
     672            dhr2(ij) = 1._wp / ( hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1) ) 
     673 
     674 
     675            ! update jl1 (removing ridged/rafted area) 
     676            a_i  (ji,jj,  jl1) = a_i  (ji,jj,  jl1) - ardg1 (ij) - arft1 (ij) 
     677            v_i  (ji,jj,  jl1) = v_i  (ji,jj,  jl1) - vrdg1 (ij) - virft (ij) 
     678            v_s  (ji,jj,  jl1) = v_s  (ji,jj,  jl1) - vsrdg (ij) - vsrft (ij) 
     679            e_s  (ji,jj,1,jl1) = e_s  (ji,jj,1,jl1) - esrdg (ij) - esrft (ij) 
     680            smv_i(ji,jj,  jl1) = smv_i(ji,jj,  jl1) - srdg1 (ij) - smrft (ij) 
     681            oa_i (ji,jj,  jl1) = oa_i (ji,jj,  jl1) - oirdg1(ij) - oirft1(ij) 
     682 
     683         END DO 
     684 
     685         !-------------------------------------------------------------------- 
     686         ! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and 
     687         !      compute ridged ice enthalpy  
     688         !-------------------------------------------------------------------- 
     689         DO jk = 1, nlay_i 
     690            DO ij = 1, icells 
     691               ji = indxi(ij) ; jj = indxj(ij) 
     692               ! heat content of ridged ice 
     693               erdg1(ij,jk) = e_i(ji,jj,jk,jl1) * afrac(ij)  
     694               eirft(ij,jk) = e_i(ji,jj,jk,jl1) * afrft(ij)                
     695                
     696               ! enthalpy of the trapped seawater (J/m2, >0) 
     697               ! clem: if sst>0, then ersw <0 (is that possible?) 
     698               ersw(ij,jk)  = - rhoic * vsw(ij) * rcp * sst_m(ji,jj) * r1_nlay_i 
     699 
     700               ! heat flux to the ocean 
     701               hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ij,jk) * r1_rdtice  ! > 0 [W.m-2] ocean->ice flux  
     702 
     703               ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean 
     704               erdg2(ij,jk) = erdg1(ij,jk) + ersw(ij,jk) 
     705 
     706               ! update jl1 
     707               e_i  (ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ij,jk) - eirft(ij,jk) 
     708 
     709            END DO 
     710         END DO 
     711 
     712         !------------------------------------------------------------------------------- 
     713         ! 4) Add area, volume, and energy of new ridge to each category jl2 
     714         !------------------------------------------------------------------------------- 
     715         DO jl2  = 1, jpl  
     716            ! over categories to which ridged/rafted ice is transferred 
     717            DO ij = 1, icells 
     718               ji = indxi(ij) ; jj = indxj(ij) 
     719 
     720               ! Compute the fraction of ridged ice area and volume going to thickness category jl2. 
     721               IF( hrmin(ji,jj,jl1) <= hi_max(jl2) .AND. hrmax(ji,jj,jl1) > hi_max(jl2-1) ) THEN 
     722                  hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) ) 
     723                  hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2)   ) 
     724                  farea    = ( hR      - hL      ) * dhr(ij)  
     725                  fvol(ij) = ( hR * hR - hL * hL ) * dhr2(ij) 
     726               ELSE 
     727                  farea    = 0._wp  
     728                  fvol(ij) = 0._wp                   
     729               ENDIF 
     730 
     731               ! Compute the fraction of rafted ice area and volume going to thickness category jl2 
     732               IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) >  hi_max(jl2-1) ) THEN 
     733                  zswitch(ij) = 1._wp 
     734               ELSE 
     735                  zswitch(ij) = 0._wp                   
     736               ENDIF 
     737 
     738               a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + ( ardg2 (ij) * farea    + arft2 (ij) * zswitch(ij) ) 
     739               oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + ( oirdg2(ij) * farea    + oirft2(ij) * zswitch(ij) ) 
     740               v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + ( vrdg2 (ij) * fvol(ij) + virft (ij) * zswitch(ij) ) 
     741               smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + ( srdg2 (ij) * fvol(ij) + smrft (ij) * zswitch(ij) ) 
     742               v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + ( vsrdg (ij) * rn_fsnowrdg * fvol(ij)  +  & 
     743                  &                                        vsrft (ij) * rn_fsnowrft * zswitch(ij) ) 
     744               e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + ( esrdg (ij) * rn_fsnowrdg * fvol(ij)  +  & 
     745                  &                                        esrft (ij) * rn_fsnowrft * zswitch(ij) ) 
     746 
     747            END DO 
     748 
     749            ! Transfer ice energy to category jl2 by ridging 
     750            DO jk = 1, nlay_i 
     751               DO ij = 1, icells 
     752                  ji = indxi(ij) ; jj = indxj(ij) 
     753                  e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + erdg2(ij,jk) * fvol(ij) + eirft(ij,jk) * zswitch(ij)                   
     754               END DO 
     755            END DO 
     756            ! 
     757         END DO ! jl2 
     758          
     759      END DO ! jl1 (deforming categories) 
     760 
     761      ! 
     762      CALL wrk_dealloc( jpij,        indxi, indxj ) 
     763      CALL wrk_dealloc( jpij,        zswitch, fvol ) 
     764      CALL wrk_dealloc( jpij,        afrac, ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 ) 
     765      CALL wrk_dealloc( jpij,        vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw, oirdg1, oirdg2 ) 
     766      CALL wrk_dealloc( jpij,        afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 ) 
     767      CALL wrk_dealloc( jpij,nlay_i, eirft, erdg1, erdg2, ersw ) 
     768      ! 
     769   END SUBROUTINE lim_itd_me_ridgeshift 
    417770 
    418771   SUBROUTINE lim_itd_me_icestrength( kstrngth ) 
    … …  
    434787      INTEGER             ::   ksmooth     ! smoothing the resistance to deformation 
    435788      INTEGER             ::   numts_rm    ! number of time steps for the P smoothing 
    436       REAL(wp)            ::   zhi, zp, z1_3  ! local scalars 
     789      REAL(wp)            ::   zp, z1_3    ! local scalars 
    437790      REAL(wp), POINTER, DIMENSION(:,:) ::   zworka   ! temporary array used here 
    438791      !!---------------------------------------------------------------------- 
    … …  
    459812               DO ji = 1, jpi 
    460813                  ! 
    461                   IF( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp ) THEN 
    462                      zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
     814                  IF( athorn(ji,jj,jl) > 0._wp ) THEN 
    463815                     !---------------------------- 
    464816                     ! PE loss from deforming ice 
    465817                     !---------------------------- 
    466                      strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * zhi * zhi 
     818                     strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * ht_i(ji,jj,jl) * ht_i(ji,jj,jl) 
    467819 
    468820                     !-------------------------- 
    469821                     ! PE gain from rafting ice 
    470822                     !-------------------------- 
    471                      strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * zhi * zhi 
     823                     strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * ht_i(ji,jj,jl) * ht_i(ji,jj,jl) 
    472824 
    473825                     !---------------------------- 
    474826                     ! PE gain from ridging ice 
    475827                     !---------------------------- 
    476                      strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) / krdg(ji,jj,jl)     & 
    477                         * z1_3 * ( hrmax(ji,jj,jl)**2 + hrmin(ji,jj,jl)**2 +  hrmax(ji,jj,jl) * hrmin(ji,jj,jl) )   
     828                     strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) * krdg(ji,jj,jl) * z1_3 *  & 
     829                        &                              ( hrmax(ji,jj,jl) * hrmax(ji,jj,jl) +         & 
     830                        &                                hrmin(ji,jj,jl) * hrmin(ji,jj,jl) +         & 
     831                        &                                hrmax(ji,jj,jl) * hrmin(ji,jj,jl) )   
    478832                        !!(a**3-b**3)/(a-b) = a*a+ab+b*b                       
    479833                  ENDIF 
    … …  
    497851         ! 
    498852      ENDIF                     ! kstrngth 
    499  
    500853      ! 
    501854      !------------------------------------------------------------------------------! 
    … …  
    503856      !------------------------------------------------------------------------------! 
    504857      ! CAN BE REMOVED 
    505       ! 
    506858      IF( ln_icestr_bvf ) THEN 
    507  
    508859         DO jj = 1, jpj 
    509860            DO ji = 1, jpi 
    … …  
    511862            END DO 
    512863         END DO 
    513  
    514864      ENDIF 
    515  
    516865      ! 
    517866      !------------------------------------------------------------------------------! 
    … …  
    558907      IF ( ksmooth == 2 ) THEN 
    559908 
    560  
    561909         CALL lbc_lnk( strength, 'T', 1. ) 
    562910 
    … …  
    565913               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp) THEN  
    566914                  numts_rm = 1 ! number of time steps for the running mean 
    567                   IF ( strp1(ji,jj) > 0.0 ) numts_rm = numts_rm + 1 
    568                   IF ( strp2(ji,jj) > 0.0 ) numts_rm = numts_rm + 1 
     915                  IF ( strp1(ji,jj) > 0._wp ) numts_rm = numts_rm + 1 
     916                  IF ( strp2(ji,jj) > 0._wp ) numts_rm = numts_rm + 1 
    569917                  zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) / numts_rm 
    570918                  strp2(ji,jj) = strp1(ji,jj) 
    … …  
    583931      ! 
    584932   END SUBROUTINE lim_itd_me_icestrength 
    585  
    586  
    587    SUBROUTINE lim_itd_me_ridgeprep 
    588       !!---------------------------------------------------------------------! 
    589       !!                ***  ROUTINE lim_itd_me_ridgeprep *** 
    590       !! 
    591       !! ** Purpose :   preparation for ridging and strength calculations 
    592       !! 
    593       !! ** Method  :   Compute the thickness distribution of the ice and open water  
    594       !!              participating in ridging and of the resulting ridges. 
    595       !!---------------------------------------------------------------------! 
    596       INTEGER ::   ji,jj, jl    ! dummy loop indices 
    597       REAL(wp) ::   Gstari, astari, zhi, hrmean, zdummy   ! local scalar 
    598       REAL(wp), POINTER, DIMENSION(:,:)   ::   zworka    ! temporary array used here 
    599       REAL(wp), POINTER, DIMENSION(:,:,:) ::   Gsum      ! Gsum(n) = sum of areas in categories 0 to n 
    600       !------------------------------------------------------------------------------! 
    601  
    602       CALL wrk_alloc( jpi,jpj, zworka ) 
    603       CALL wrk_alloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 ) 
    604  
    605       Gstari     = 1.0/rn_gstar     
    606       astari     = 1.0/rn_astar     
    607       aksum(:,:)    = 0.0 
    608       athorn(:,:,:) = 0.0 
    609       aridge(:,:,:) = 0.0 
    610       araft (:,:,:) = 0.0 
    611       hrmin(:,:,:)  = 0.0  
    612       hrmax(:,:,:)  = 0.0  
    613       hraft(:,:,:)  = 0.0  
    614       krdg (:,:,:)  = 1.0 
    615  
    616       !     ! Zero out categories with very small areas 
    617       CALL lim_var_zapsmall 
    618  
    619       !------------------------------------------------------------------------------! 
    620       ! 1) Participation function  
    621       !------------------------------------------------------------------------------! 
    622  
    623       ! Compute total area of ice plus open water. 
    624       ! This is in general not equal to one because of divergence during transport 
    625       asum(:,:) = ato_i(:,:) 
    626       DO jl = 1, jpl 
    627          asum(:,:) = asum(:,:) + a_i(:,:,jl) 
    628       END DO 
    629  
    630       ! Compute cumulative thickness distribution function 
    631       ! Compute the cumulative thickness distribution function Gsum, 
    632       ! where Gsum(n) is the fractional area in categories 0 to n. 
    633       ! initial value (in h = 0) equals open water area 
    634  
    635       Gsum(:,:,-1) = 0._wp 
    636       Gsum(:,:,0 ) = ato_i(:,:) 
    637  
    638       ! for each value of h, you have to add ice concentration then 
    639       DO jl = 1, jpl 
    640          Gsum(:,:,jl) = Gsum(:,:,jl-1) + a_i(:,:,jl) 
    641       END DO 
    642  
    643       ! Normalize the cumulative distribution to 1 
    644       zworka(:,:) = 1._wp / Gsum(:,:,jpl) 
    645       DO jl = 0, jpl 
    646          Gsum(:,:,jl) = Gsum(:,:,jl) * zworka(:,:) 
    647       END DO 
    648  
    649       ! 1.3 Compute participation function a(h) = b(h).g(h) (athorn) 
    650       !-------------------------------------------------------------------------------------------------- 
    651       ! Compute the participation function athorn; this is analogous to 
    652       ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975). 
    653       ! area lost from category n due to ridging/closing 
    654       ! athorn(n)   = total area lost due to ridging/closing 
    655       ! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).  
    656       ! 
    657       ! The expressions for athorn are found by integrating b(h)g(h) between 
    658       ! the category boundaries. 
    659       !----------------------------------------------------------------- 
    660  
    661       IF( nn_partfun == 0 ) THEN       !--- Linear formulation (Thorndike et al., 1975) 
    662          DO jl = 0, jpl     
    663             DO jj = 1, jpj  
    664                DO ji = 1, jpi 
    665                   IF( Gsum(ji,jj,jl) < rn_gstar) THEN 
    666                      athorn(ji,jj,jl) = Gstari * ( Gsum(ji,jj,jl) - Gsum(ji,jj,jl-1) ) * & 
    667                         &                        ( 2.0 - (Gsum(ji,jj,jl-1) + Gsum(ji,jj,jl) ) * Gstari ) 
    668                   ELSEIF (Gsum(ji,jj,jl-1) < rn_gstar) THEN 
    669                      athorn(ji,jj,jl) = Gstari * ( rn_gstar - Gsum(ji,jj,jl-1) ) *  & 
    670                         &                        ( 2.0 - ( Gsum(ji,jj,jl-1) + rn_gstar ) * Gstari ) 
    671                   ELSE 
    672                      athorn(ji,jj,jl) = 0.0 
    673                   ENDIF 
    674                END DO 
    675             END DO 
    676          END DO 
    677  
    678       ELSE                             !--- Exponential, more stable formulation (Lipscomb et al, 2007) 
    679          !                         
    680          zdummy = 1._wp / ( 1._wp - EXP(-astari) )        ! precompute exponential terms using Gsum as a work array 
    681          DO jl = -1, jpl 
    682             Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy 
    683          END DO 
    684          DO jl = 0, jpl 
    685              athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl) 
    686          END DO 
    687          ! 
    688       ENDIF 
    689  
    690       IF( ln_rafting ) THEN      ! Ridging and rafting ice participation functions 
    691          ! 
    692          DO jl = 1, jpl 
    693             DO jj = 1, jpj  
    694                DO ji = 1, jpi 
    695                   IF ( athorn(ji,jj,jl) > 0._wp ) THEN 
    696 !!gm  TANH( -X ) = - TANH( X )  so can be computed only 1 time.... 
    697                      aridge(ji,jj,jl) = ( TANH (  rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl) 
    698                      araft (ji,jj,jl) = ( TANH ( -rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl) 
    699                      IF ( araft(ji,jj,jl) < epsi06 )   araft(ji,jj,jl)  = 0._wp 
    700                      aridge(ji,jj,jl) = MAX( athorn(ji,jj,jl) - araft(ji,jj,jl), 0.0 ) 
    701                   ENDIF 
    702                END DO 
    703             END DO 
    704          END DO 
    705  
    706       ELSE 
    707          ! 
    708          DO jl = 1, jpl 
    709             aridge(:,:,jl) = athorn(:,:,jl) 
    710          END DO 
    711          ! 
    712       ENDIF 
    713  
    714       IF( ln_rafting ) THEN 
    715  
    716          IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) > epsi10 .AND. lwp ) THEN 
    717             DO jl = 1, jpl 
    718                DO jj = 1, jpj 
    719                   DO ji = 1, jpi 
    720                      IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) > epsi10 ) THEN 
    721                         WRITE(numout,*) ' ALERTE 96 : wrong participation function ... ' 
    722                         WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl 
    723                         WRITE(numout,*) ' lat, lon   : ', gphit(ji,jj), glamt(ji,jj) 
    724                         WRITE(numout,*) ' aridge     : ', aridge(ji,jj,1:jpl) 
    725                         WRITE(numout,*) ' araft      : ', araft(ji,jj,1:jpl) 
    726                         WRITE(numout,*) ' athorn     : ', athorn(ji,jj,1:jpl) 
    727                      ENDIF 
    728                   END DO 
    729                END DO 
    730             END DO 
    731          ENDIF 
    732  
    733       ENDIF 
    734  
    735       !----------------------------------------------------------------- 
    736       ! 2) Transfer function 
    737       !----------------------------------------------------------------- 
    738       ! Compute max and min ridged ice thickness for each ridging category. 
    739       ! Assume ridged ice is uniformly distributed between hrmin and hrmax. 
    740       !  
    741       ! This parameterization is a modified version of Hibler (1980). 
    742       ! The mean ridging thickness, hrmean, is proportional to hi^(0.5) 
    743       !  and for very thick ridging ice must be >= krdgmin*hi 
    744       ! 
    745       ! The minimum ridging thickness, hrmin, is equal to 2*hi  
    746       !  (i.e., rafting) and for very thick ridging ice is 
    747       !  constrained by hrmin <= (hrmean + hi)/2. 
    748       !  
    749       ! The maximum ridging thickness, hrmax, is determined by 
    750       !  hrmean and hrmin. 
    751       ! 
    752       ! These modifications have the effect of reducing the ice strength 
    753       ! (relative to the Hibler formulation) when very thick ice is 
    754       ! ridging. 
    755       ! 
    756       ! aksum = net area removed/ total area removed 
    757       ! where total area removed = area of ice that ridges 
    758       !         net area removed = total area removed - area of new ridges 
    759       !----------------------------------------------------------------- 
    760  
    761       ! Transfer function 
    762       DO jl = 1, jpl !all categories have a specific transfer function 
    763          DO jj = 1, jpj 
    764             DO ji = 1, jpi 
    765  
    766                IF (a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0.0 ) THEN 
    767                   zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
    768                   hrmean          = MAX(SQRT(rn_hstar*zhi), zhi*krdgmin) 
    769                   hrmin(ji,jj,jl) = MIN(2.0*zhi, 0.5*(hrmean + zhi)) 
    770                   hrmax(ji,jj,jl) = 2.0*hrmean - hrmin(ji,jj,jl) 
    771                   hraft(ji,jj,jl) = kraft*zhi 
    772                   krdg(ji,jj,jl)  = hrmean / zhi 
    773                ELSE 
    774                   hraft(ji,jj,jl) = 0.0 
    775                   hrmin(ji,jj,jl) = 0.0  
    776                   hrmax(ji,jj,jl) = 0.0  
    777                   krdg (ji,jj,jl) = 1.0 
    778                ENDIF 
    779  
    780             END DO 
    781          END DO 
    782       END DO 
    783  
    784       ! Normalization factor : aksum, ensures mass conservation 
    785       aksum(:,:) = athorn(:,:,0) 
    786       DO jl = 1, jpl  
    787          aksum(:,:)    = aksum(:,:) + aridge(:,:,jl) * ( 1._wp - 1._wp / krdg(:,:,jl) )    & 
    788             &                       + araft (:,:,jl) * ( 1._wp - 1._wp / kraft        ) 
    789       END DO 
    790       ! 
    791       CALL wrk_dealloc( jpi,jpj, zworka ) 
    792       CALL wrk_dealloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 ) 
    793       ! 
    794    END SUBROUTINE lim_itd_me_ridgeprep 
    795  
    796  
    797    SUBROUTINE lim_itd_me_ridgeshift( opning, closing_gross, msnow_mlt, esnow_mlt ) 
    798       !!---------------------------------------------------------------------- 
    799       !!                ***  ROUTINE lim_itd_me_icestrength *** 
    800       !! 
    801       !! ** Purpose :   shift ridging ice among thickness categories of ice thickness 
    802       !! 
    803       !! ** Method  :   Remove area, volume, and energy from each ridging category 
    804       !!              and add to thicker ice categories. 
    805       !!---------------------------------------------------------------------- 
    806       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   opning         ! rate of opening due to divergence/shear 
    807       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   closing_gross  ! rate at which area removed, excluding area of new ridges 
    808       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   msnow_mlt      ! mass of snow added to ocean (kg m-2) 
    809       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   esnow_mlt      ! energy needed to melt snow in ocean (J m-2) 
    810       ! 
    811       CHARACTER (len=80) ::   fieldid   ! field identifier 
    812       LOGICAL, PARAMETER ::   l_conservation_check = .true.  ! if true, check conservation (useful for debugging) 
    813       ! 
    814       INTEGER ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices 
    815       INTEGER ::   ij                ! horizontal index, combines i and j loops 
    816       INTEGER ::   icells            ! number of cells with aicen > puny 
    817       REAL(wp) ::   hL, hR, farea, ztmelts    ! left and right limits of integration 
    818  
    819       INTEGER , POINTER, DIMENSION(:) ::   indxi, indxj   ! compressed indices 
    820  
    821       REAL(wp), POINTER, DIMENSION(:,:) ::   vice_init, vice_final   ! ice volume summed over categories 
    822       REAL(wp), POINTER, DIMENSION(:,:) ::   eice_init, eice_final   ! ice energy summed over layers 
    823  
    824       REAL(wp), POINTER, DIMENSION(:,:,:) ::   aicen_init, vicen_init   ! ice  area    & volume before ridging 
    825       REAL(wp), POINTER, DIMENSION(:,:,:) ::   vsnwn_init, esnwn_init   ! snow volume  & energy before ridging 
    826       REAL(wp), POINTER, DIMENSION(:,:,:) ::   smv_i_init, oa_i_init    ! ice salinity & age    before ridging 
    827  
    828       REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   eicen_init        ! ice energy before ridging 
    829  
    830       REAL(wp), POINTER, DIMENSION(:,:) ::   afrac , fvol     ! fraction of category area ridged & new ridge volume going to n2 
    831       REAL(wp), POINTER, DIMENSION(:,:) ::   ardg1 , ardg2    ! area of ice ridged & new ridges 
    832       REAL(wp), POINTER, DIMENSION(:,:) ::   vsrdg , esrdg    ! snow volume & energy of ridging ice 
    833       REAL(wp), POINTER, DIMENSION(:,:) ::   dhr   , dhr2     ! hrmax - hrmin  &  hrmax^2 - hrmin^2 
    834  
    835       REAL(wp), POINTER, DIMENSION(:,:) ::   vrdg1   ! volume of ice ridged 
    836       REAL(wp), POINTER, DIMENSION(:,:) ::   vrdg2   ! volume of new ridges 
    837       REAL(wp), POINTER, DIMENSION(:,:) ::   vsw     ! volume of seawater trapped into ridges 
    838       REAL(wp), POINTER, DIMENSION(:,:) ::   srdg1   ! sal*volume of ice ridged 
    839       REAL(wp), POINTER, DIMENSION(:,:) ::   srdg2   ! sal*volume of new ridges 
    840       REAL(wp), POINTER, DIMENSION(:,:) ::   smsw    ! sal*volume of water trapped into ridges 
    841       REAL(wp), POINTER, DIMENSION(:,:) ::   oirdg1, oirdg2   ! ice age of ice ridged 
    842  
    843       REAL(wp), POINTER, DIMENSION(:,:) ::   afrft            ! fraction of category area rafted 
    844       REAL(wp), POINTER, DIMENSION(:,:) ::   arft1 , arft2    ! area of ice rafted and new rafted zone 
    845       REAL(wp), POINTER, DIMENSION(:,:) ::   virft , vsrft    ! ice & snow volume of rafting ice 
    846       REAL(wp), POINTER, DIMENSION(:,:) ::   esrft , smrft    ! snow energy & salinity of rafting ice 
    847       REAL(wp), POINTER, DIMENSION(:,:) ::   oirft1, oirft2   ! ice age of ice rafted 
    848  
    849       REAL(wp), POINTER, DIMENSION(:,:,:) ::   eirft      ! ice energy of rafting ice 
    850       REAL(wp), POINTER, DIMENSION(:,:,:) ::   erdg1      ! enth*volume of ice ridged 
    851       REAL(wp), POINTER, DIMENSION(:,:,:) ::   erdg2      ! enth*volume of new ridges 
    852       REAL(wp), POINTER, DIMENSION(:,:,:) ::   ersw       ! enth of water trapped into ridges 
    853       !!---------------------------------------------------------------------- 
    854  
    855       CALL wrk_alloc( (jpi+1)*(jpj+1),       indxi, indxj ) 
    856       CALL wrk_alloc( jpi, jpj,              vice_init, vice_final, eice_init, eice_final ) 
    857       CALL wrk_alloc( jpi, jpj,              afrac, fvol , ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 ) 
    858       CALL wrk_alloc( jpi, jpj,              vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw, oirdg1, oirdg2 ) 
    859       CALL wrk_alloc( jpi, jpj,              afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 ) 
    860       CALL wrk_alloc( jpi, jpj, jpl,         aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init ) 
    861       CALL wrk_alloc( jpi, jpj, nlay_i,      eirft, erdg1, erdg2, ersw ) 
    862       CALL wrk_alloc( jpi, jpj, nlay_i, jpl, eicen_init ) 
    863  
    864       ! Conservation check 
    865       eice_init(:,:) = 0._wp 
    866  
    867       IF( con_i ) THEN 
    868          CALL lim_column_sum        (jpl,    v_i,       vice_init ) 
    869          CALL lim_column_sum_energy (jpl, nlay_i,  e_i, eice_init ) 
    870          DO ji = mi0(iiceprt), mi1(iiceprt) 
    871             DO jj = mj0(jiceprt), mj1(jiceprt) 
    872                WRITE(numout,*) ' vice_init  : ', vice_init(ji,jj) 
    873                WRITE(numout,*) ' eice_init  : ', eice_init(ji,jj) 
    874             END DO 
    875          END DO 
    876       ENDIF 
    877  
    878       !------------------------------------------------------------------------------- 
    879       ! 1) Compute change in open water area due to closing and opening. 
    880       !------------------------------------------------------------------------------- 
    881       DO jj = 1, jpj 
    882          DO ji = 1, jpi 
    883             ato_i(ji,jj) = ato_i(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice        & 
    884                &                        + opning(ji,jj)                          * rdt_ice 
    885             IF    ( ato_i(ji,jj) < -epsi10 ) THEN    ! there is a bug 
    886                IF(lwp)   WRITE(numout,*) 'Ridging error: ato_i < 0 -- ato_i : ',ato_i(ji,jj) 
    887             ELSEIF( ato_i(ji,jj) < 0._wp   ) THEN    ! roundoff error 
    888                ato_i(ji,jj) = 0._wp 
    889             ENDIF 
    890          END DO 
    891       END DO 
    892  
    893       !----------------------------------------------------------------- 
    894       ! 2) Save initial state variables 
    895       !----------------------------------------------------------------- 
    896       aicen_init(:,:,:)   = a_i  (:,:,:) 
    897       vicen_init(:,:,:)   = v_i  (:,:,:) 
    898       vsnwn_init(:,:,:)   = v_s  (:,:,:) 
    899       smv_i_init(:,:,:)   = smv_i(:,:,:) 
    900       esnwn_init(:,:,:)   = e_s  (:,:,1,:) 
    901       eicen_init(:,:,:,:) = e_i  (:,:,:,:) 
    902       oa_i_init (:,:,:)   = oa_i (:,:,:) 
    903  
    904       ! 
    905       !----------------------------------------------------------------- 
    906       ! 3) Pump everything from ice which is being ridged / rafted 
    907       !----------------------------------------------------------------- 
    908       ! Compute the area, volume, and energy of ice ridging in each 
    909       ! category, along with the area of the resulting ridge. 
    910  
    911       DO jl1 = 1, jpl !jl1 describes the ridging category 
    912  
    913          !------------------------------------------------ 
    914          ! 3.1) Identify grid cells with nonzero ridging 
    915          !------------------------------------------------ 
    916  
    917          icells = 0 
    918          DO jj = 1, jpj 
    919             DO ji = 1, jpi 
    920                IF( aicen_init(ji,jj,jl1) > epsi10 .AND. athorn(ji,jj,jl1) > 0._wp  & 
    921                   &   .AND. closing_gross(ji,jj) > 0._wp ) THEN 
    922                   icells = icells + 1 
    923                   indxi(icells) = ji 
    924                   indxj(icells) = jj 
    925                ENDIF 
    926             END DO 
    927          END DO 
    928  
    929          DO ij = 1, icells 
    930             ji = indxi(ij) 
    931             jj = indxj(ij) 
    932  
    933             !-------------------------------------------------------------------- 
    934             ! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2) 
    935             !-------------------------------------------------------------------- 
    936  
    937             ardg1(ji,jj) = aridge(ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice 
    938             arft1(ji,jj) = araft (ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice 
    939             ardg2(ji,jj) = ardg1(ji,jj) / krdg(ji,jj,jl1) 
    940             arft2(ji,jj) = arft1(ji,jj) / kraft 
    941  
    942             !--------------------------------------------------------------- 
    943             ! 3.3) Compute ridging /rafting fractions, make sure afrac <=1  
    944             !--------------------------------------------------------------- 
    945  
    946             afrac(ji,jj) = ardg1(ji,jj) / aicen_init(ji,jj,jl1) !ridging 
    947             afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting 
    948  
    949             IF( afrac(ji,jj) > kamax + epsi10 ) THEN  ! there is a bug 
    950                IF(lwp)   WRITE(numout,*) ' ardg > a_i -- ardg, aicen_init : ', ardg1(ji,jj), aicen_init(ji,jj,jl1) 
    951             ELSEIF( afrac(ji,jj) > kamax ) THEN       ! roundoff error 
    952                afrac(ji,jj) = kamax 
    953             ENDIF 
    954  
    955             IF( afrft(ji,jj) > kamax + epsi10 ) THEN ! there is a bug 
    956                IF(lwp)   WRITE(numout,*) ' arft > a_i -- arft, aicen_init : ', arft1(ji,jj), aicen_init(ji,jj,jl1)  
    957             ELSEIF( afrft(ji,jj) > kamax) THEN       ! roundoff error 
    958                afrft(ji,jj) = kamax 
    959             ENDIF 
    960  
    961             !-------------------------------------------------------------------------- 
    962             ! 3.4) Subtract area, volume, and energy from ridging  
    963             !     / rafting category n1. 
    964             !-------------------------------------------------------------------------- 
    965             vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) 
    966             vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + rn_por_rdg ) 
    967             vsw  (ji,jj) = vrdg1(ji,jj) * rn_por_rdg 
    968  
    969             vsrdg (ji,jj) = vsnwn_init(ji,jj,jl1) * afrac(ji,jj) 
    970             esrdg (ji,jj) = esnwn_init(ji,jj,jl1) * afrac(ji,jj) 
    971             srdg1 (ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) 
    972             oirdg1(ji,jj) = oa_i_init (ji,jj,jl1) * afrac(ji,jj) 
    973             oirdg2(ji,jj) = oa_i_init (ji,jj,jl1) * afrac(ji,jj) / krdg(ji,jj,jl1)  
    974  
    975             ! rafting volumes, heat contents ... 
    976             virft (ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj) 
    977             vsrft (ji,jj) = vsnwn_init(ji,jj,jl1) * afrft(ji,jj) 
    978             esrft (ji,jj) = esnwn_init(ji,jj,jl1) * afrft(ji,jj) 
    979             smrft (ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj)  
    980             oirft1(ji,jj) = oa_i_init (ji,jj,jl1) * afrft(ji,jj)  
    981             oirft2(ji,jj) = oa_i_init (ji,jj,jl1) * afrft(ji,jj) / kraft  
    982  
    983             ! substract everything 
    984             a_i(ji,jj,jl1)   = a_i(ji,jj,jl1)   - ardg1 (ji,jj) - arft1 (ji,jj) 
    985             v_i(ji,jj,jl1)   = v_i(ji,jj,jl1)   - vrdg1 (ji,jj) - virft (ji,jj) 
    986             v_s(ji,jj,jl1)   = v_s(ji,jj,jl1)   - vsrdg (ji,jj) - vsrft (ji,jj) 
    987             e_s(ji,jj,1,jl1) = e_s(ji,jj,1,jl1) - esrdg (ji,jj) - esrft (ji,jj) 
    988             smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - srdg1 (ji,jj) - smrft (ji,jj) 
    989             oa_i(ji,jj,jl1)  = oa_i(ji,jj,jl1)  - oirdg1(ji,jj) - oirft1(ji,jj) 
    990  
    991             !----------------------------------------------------------------- 
    992             ! 3.5) Compute properties of new ridges 
    993             !----------------------------------------------------------------- 
    994             !--------- 
    995             ! Salinity 
    996             !--------- 
    997             smsw(ji,jj)  = vsw(ji,jj) * sss_m(ji,jj)                      ! salt content of seawater frozen in voids !! MV HC2014 
    998             srdg2(ji,jj) = srdg1(ji,jj) + smsw(ji,jj)                     ! salt content of new ridge 
    999  
    1000             !srdg2(ji,jj) = MIN( rn_simax * vrdg2(ji,jj) , zsrdg2 )         ! impose a maximum salinity 
    1001              
    1002             sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ji,jj) * rhoic * r1_rdtice 
    1003             wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ji,jj) * rhoic * r1_rdtice   ! increase in ice volume du to seawater frozen in voids              
    1004  
    1005             !------------------------------------             
    1006             ! 3.6 Increment ridging diagnostics 
    1007             !------------------------------------             
    1008  
    1009             !        jl1 looping 1-jpl 
    1010             !           ij looping 1-icells 
    1011  
    1012             dardg1dt   (ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj) 
    1013             dardg2dt   (ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj) 
    1014             opening    (ji,jj) = opening (ji,jj) + opning(ji,jj) * rdt_ice 
    1015  
    1016             IF( con_i )   vice_init(ji,jj) = vice_init(ji,jj) + vrdg2(ji,jj) - vrdg1(ji,jj) 
    1017  
    1018             !------------------------------------------             
    1019             ! 3.7 Put the snow somewhere in the ocean 
    1020             !------------------------------------------             
    1021             !  Place part of the snow lost by ridging into the ocean.  
    1022             !  Note that esnow_mlt < 0; the ocean must cool to melt snow. 
    1023             !  If the ocean temp = Tf already, new ice must grow. 
    1024             !  During the next time step, thermo_rates will determine whether 
    1025             !  the ocean cools or new ice grows. 
    1026             !        jl1 looping 1-jpl 
    1027             !           ij looping 1-icells 
    1028  
    1029             msnow_mlt(ji,jj) = msnow_mlt(ji,jj) + rhosn*vsrdg(ji,jj)*(1.0-rn_fsnowrdg)   &   ! rafting included 
    1030                &                                + rhosn*vsrft(ji,jj)*(1.0-rn_fsnowrft) 
    1031  
    1032             ! in J/m2 (same as e_s) 
    1033             esnow_mlt(ji,jj) = esnow_mlt(ji,jj) - esrdg(ji,jj)*(1.0-rn_fsnowrdg)         &   !rafting included 
    1034                &                                - esrft(ji,jj)*(1.0-rn_fsnowrft)           
    1035  
    1036             !----------------------------------------------------------------- 
    1037             ! 3.8 Compute quantities used to apportion ice among categories 
    1038             ! in the n2 loop below 
    1039             !----------------------------------------------------------------- 
    1040  
    1041             !        jl1 looping 1-jpl 
    1042             !           ij looping 1-icells 
    1043  
    1044             dhr (ji,jj) = hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) 
    1045             dhr2(ji,jj) = hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1) 
    1046  
    1047          END DO 
    1048  
    1049          !-------------------------------------------------------------------- 
    1050          ! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and 
    1051          !      compute ridged ice enthalpy  
    1052          !-------------------------------------------------------------------- 
    1053          DO jk = 1, nlay_i 
    1054             DO ij = 1, icells 
    1055                ji = indxi(ij) 
    1056                jj = indxj(ij) 
    1057                ! heat content of ridged ice 
    1058                erdg1(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj)  
    1059                eirft(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrft(ji,jj) 
    1060                e_i  (ji,jj,jk,jl1)  = e_i(ji,jj,jk,jl1) - erdg1(ji,jj,jk) - eirft(ji,jj,jk) 
    1061                 
    1062                 
    1063                ! enthalpy of the trapped seawater (J/m2, >0) 
    1064                ! clem: if sst>0, then ersw <0 (is that possible?) 
    1065                ersw(ji,jj,jk)   = - rhoic * vsw(ji,jj) * rcp * sst_m(ji,jj) * r1_nlay_i 
    1066  
    1067                ! heat flux to the ocean 
    1068                hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ji,jj,jk) * r1_rdtice  ! > 0 [W.m-2] ocean->ice flux  
    1069  
    1070                ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean 
    1071                erdg2(ji,jj,jk)  = erdg1(ji,jj,jk) + ersw(ji,jj,jk) 
    1072  
    1073             END DO 
    1074          END DO 
    1075  
    1076  
    1077          IF( con_i ) THEN 
    1078             DO jk = 1, nlay_i 
    1079                DO ij = 1, icells 
    1080                   ji = indxi(ij) 
    1081                   jj = indxj(ij) 
    1082                   eice_init(ji,jj) = eice_init(ji,jj) + erdg2(ji,jj,jk) - erdg1(ji,jj,jk) 
    1083                END DO 
    1084             END DO 
    1085          ENDIF 
    1086  
    1087          !------------------------------------------------------------------------------- 
    1088          ! 4) Add area, volume, and energy of new ridge to each category jl2 
    1089          !------------------------------------------------------------------------------- 
    1090          !        jl1 looping 1-jpl 
    1091          DO jl2  = 1, jpl  
    1092             ! over categories to which ridged ice is transferred 
    1093             DO ij = 1, icells 
    1094                ji = indxi(ij) 
    1095                jj = indxj(ij) 
    1096  
    1097                ! Compute the fraction of ridged ice area and volume going to  
    1098                ! thickness category jl2. 
    1099                ! Transfer area, volume, and energy accordingly. 
    1100  
    1101                IF( hrmin(ji,jj,jl1) >= hi_max(jl2) .OR. hrmax(ji,jj,jl1) <= hi_max(jl2-1) ) THEN 
    1102                   hL = 0._wp 
    1103                   hR = 0._wp 
    1104                ELSE 
    1105                   hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) ) 
    1106                   hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2)   ) 
    1107                ENDIF 
    1108  
    1109                ! fraction of ridged ice area and volume going to n2 
    1110                farea = ( hR - hL ) / dhr(ji,jj)  
    1111                fvol(ji,jj) = ( hR*hR - hL*hL ) / dhr2(ji,jj) 
    1112  
    1113                a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + ardg2 (ji,jj) * farea 
    1114                v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + vrdg2 (ji,jj) * fvol(ji,jj) 
    1115                v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + vsrdg (ji,jj) * fvol(ji,jj) * rn_fsnowrdg 
    1116                e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + esrdg (ji,jj) * fvol(ji,jj) * rn_fsnowrdg 
    1117                smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + srdg2 (ji,jj) * fvol(ji,jj) 
    1118                oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirdg2(ji,jj) * farea 
    1119  
    1120             END DO 
    1121  
    1122             ! Transfer ice energy to category jl2 by ridging 
    1123             DO jk = 1, nlay_i 
    1124                DO ij = 1, icells 
    1125                   ji = indxi(ij) 
    1126                   jj = indxj(ij) 
    1127                   e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + fvol(ji,jj) * erdg2(ji,jj,jk) 
    1128                END DO 
    1129             END DO 
    1130             ! 
    1131          END DO                 ! jl2 (new ridges)             
    1132  
    1133          DO jl2 = 1, jpl  
    1134  
    1135             DO ij = 1, icells 
    1136                ji = indxi(ij) 
    1137                jj = indxj(ij) 
    1138                ! Compute the fraction of rafted ice area and volume going to  
    1139                ! thickness category jl2, transfer area, volume, and energy accordingly. 
    1140                ! 
    1141                IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) >  hi_max(jl2-1) ) THEN 
    1142                   a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + arft2 (ji,jj) 
    1143                   v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + virft (ji,jj) 
    1144                   v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + vsrft (ji,jj) * rn_fsnowrft 
    1145                   e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + esrft (ji,jj) * rn_fsnowrft 
    1146                   smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + smrft (ji,jj)     
    1147                   oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirft2(ji,jj) 
    1148                ENDIF 
    1149                ! 
    1150             END DO 
    1151  
    1152             ! Transfer rafted ice energy to category jl2  
    1153             DO jk = 1, nlay_i 
    1154                DO ij = 1, icells 
    1155                   ji = indxi(ij) 
    1156                   jj = indxj(ij) 
    1157                   IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) > hi_max(jl2-1)  ) THEN 
    1158                      e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + eirft(ji,jj,jk) 
    1159                   ENDIF 
    1160                END DO 
    1161             END DO 
    1162  
    1163          END DO 
    1164  
    1165       END DO ! jl1 (deforming categories) 
    1166  
    1167       ! Conservation check 
    1168       IF ( con_i ) THEN 
    1169          CALL lim_column_sum (jpl,   v_i, vice_final) 
    1170          fieldid = ' v_i : limitd_me ' 
    1171          CALL lim_cons_check (vice_init, vice_final, 1.0e-6, fieldid)  
    1172  
    1173          CALL lim_column_sum_energy (jpl, nlay_i,  e_i, eice_final ) 
    1174          fieldid = ' e_i : limitd_me ' 
    1175          CALL lim_cons_check (eice_init, eice_final, 1.0e-2, fieldid)  
    1176  
    1177          DO ji = mi0(iiceprt), mi1(iiceprt) 
    1178             DO jj = mj0(jiceprt), mj1(jiceprt) 
    1179                WRITE(numout,*) ' vice_init  : ', vice_init (ji,jj) 
    1180                WRITE(numout,*) ' vice_final : ', vice_final(ji,jj) 
    1181                WRITE(numout,*) ' eice_init  : ', eice_init (ji,jj) 
    1182                WRITE(numout,*) ' eice_final : ', eice_final(ji,jj) 
    1183             END DO 
    1184          END DO 
    1185       ENDIF 
    1186       ! 
    1187       CALL wrk_dealloc( (jpi+1)*(jpj+1),        indxi, indxj ) 
    1188       CALL wrk_dealloc( jpi, jpj,               vice_init, vice_final, eice_init, eice_final ) 
    1189       CALL wrk_dealloc( jpi, jpj,               afrac, fvol , ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 ) 
    1190       CALL wrk_dealloc( jpi, jpj,               vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw, oirdg1, oirdg2 ) 
    1191       CALL wrk_dealloc( jpi, jpj,               afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 ) 
    1192       CALL wrk_dealloc( jpi, jpj, jpl,          aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init ) 
    1193       CALL wrk_dealloc( jpi, jpj, nlay_i,       eirft, erdg1, erdg2, ersw ) 
    1194       CALL wrk_dealloc( jpi, jpj, nlay_i, jpl,  eicen_init ) 
    1195       ! 
    1196    END SUBROUTINE lim_itd_me_ridgeshift 
    1197933 
    1198934   SUBROUTINE lim_itd_me_init 

  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90

    r6333 r6436  
    125125      IF( iom_use('emp_ice' ) )  CALL iom_put( "emp_ice"  , emp_ice(:,:) )   ! emp over ice   (taking into account the snow blown away from the ice) 
    126126 
    127       ! clem 2016: albedo output 
     127      ! albedo output 
    128128      CALL wrk_alloc( jpi,jpj, zalb )     
    129129 
    … …  
    158158            hfx_out(ji,jj) = hfx_out(ji,jj) + zqmass + zqsr 
    159159 
    160             ! Add the residual from heat diffusion equation (W.m-2) 
    161             !------------------------------------------------------- 
    162             hfx_out(ji,jj) = hfx_out(ji,jj) + hfx_err_dif(ji,jj) 
     160            ! Add the residual from heat diffusion equation and sublimation (W.m-2) 
     161            !---------------------------------------------------------------------- 
     162            hfx_out(ji,jj) = hfx_out(ji,jj) + hfx_err_dif(ji,jj) +   & 
     163               &           ( hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) ) 
    163164 
    164165            ! New qsr and qns used to compute the oceanic heat flux at the next time step 
    165             !--------------------------------------------------- 
     166            !---------------------------------------------------------------------------- 
    166167            qsr(ji,jj) = zqsr                                       
    167168            qns(ji,jj) = hfx_out(ji,jj) - zqsr               
    … …  
    183184 
    184185            ! mass flux at the ocean/ice interface 
    185             fmmflx(ji,jj) = - ( wfx_ice(ji,jj) + wfx_snw(ji,jj) ) * r1_rdtice  ! F/M mass flux save at least for biogeochemical model 
    186             emp(ji,jj)    = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj)   ! mass flux + F/M mass flux (always ice/ocean mass exchange) 
    187              
     186            fmmflx(ji,jj) = - ( wfx_ice(ji,jj) + wfx_snw(ji,jj) + wfx_err_sub(ji,jj) )              ! F/M mass flux save at least for biogeochemical model 
     187            emp(ji,jj)    = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj)   ! mass flux + F/M mass flux (always ice/ocean mass exchange) 
    188188         END DO 
    189189      END DO 
    … …  
    193193      !------------------------------------------! 
    194194      sfx(:,:) = sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:)   & 
    195          &     + sfx_res(:,:) + sfx_dyn(:,:) + sfx_bri(:,:) 
     195         &     + sfx_res(:,:) + sfx_dyn(:,:) + sfx_bri(:,:) + sfx_sub(:,:) 
    196196 
    197197      !-------------------------------------------------------------! 

  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r5781 r6436  
    461461 
    462462      DO ji = kideb, kiut 
    463          zdh_mel = MIN( 0._wp, dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) ) 
     463         zdh_mel = MIN( 0._wp, dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) + dh_i_sub(ji) ) 
    464464         IF( zdh_mel < 0._wp .AND. a_i_1d(ji) > 0._wp )  THEN 
    465465            zvi          = a_i_1d(ji) * ht_i_1d(ji) 
    … …  
    470470            zda_mel     = rswitch * a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi20 ) ) 
    471471            a_i_1d(ji)  = MAX( epsi20, a_i_1d(ji) + zda_mel )  
    472              ! adjust thickness 
     472            ! adjust thickness 
    473473            ht_i_1d(ji) = zvi / a_i_1d(ji)             
    474474            ht_s_1d(ji) = zvs / a_i_1d(ji)             
    … …  
    514514          
    515515         CALL tab_2d_1d( nbpb, qprec_ice_1d(1:nbpb), qprec_ice(:,:) , jpi, jpj, npb(1:nbpb) ) 
     516         CALL tab_2d_1d( nbpb, qevap_ice_1d(1:nbpb), qevap_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 
    516517         CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb), qsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 
    517518         CALL tab_2d_1d( nbpb, fr1_i0_1d  (1:nbpb), fr1_i0          , jpi, jpj, npb(1:nbpb) ) 
    … …  
    543544         CALL tab_2d_1d( nbpb, sfx_bri_1d (1:nbpb), sfx_bri         , jpi, jpj, npb(1:nbpb) ) 
    544545         CALL tab_2d_1d( nbpb, sfx_res_1d (1:nbpb), sfx_res         , jpi, jpj, npb(1:nbpb) ) 
    545           
     546         CALL tab_2d_1d( nbpb, sfx_sub_1d (1:nbpb), sfx_sub         , jpi, jpj,npb(1:nbpb) ) 
     547  
    546548         CALL tab_2d_1d( nbpb, hfx_thd_1d (1:nbpb), hfx_thd         , jpi, jpj, npb(1:nbpb) ) 
    547549         CALL tab_2d_1d( nbpb, hfx_spr_1d (1:nbpb), hfx_spr         , jpi, jpj, npb(1:nbpb) ) 
    … …  
    593595         CALL tab_1d_2d( nbpb, sfx_res       , npb, sfx_res_1d(1:nbpb)   , jpi, jpj ) 
    594596         CALL tab_1d_2d( nbpb, sfx_bri       , npb, sfx_bri_1d(1:nbpb)   , jpi, jpj ) 
    595           
     597         CALL tab_1d_2d( nbpb, sfx_sub       , npb, sfx_sub_1d(1:nbpb)   , jpi, jpj )         
     598  
    596599         CALL tab_1d_2d( nbpb, hfx_thd       , npb, hfx_thd_1d(1:nbpb)   , jpi, jpj ) 
    597600         CALL tab_1d_2d( nbpb, hfx_spr       , npb, hfx_spr_1d(1:nbpb)   , jpi, jpj ) 

  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

    r5781 r6436  
    7474 
    7575      REAL(wp) ::   ztmelts             ! local scalar 
    76       REAL(wp) ::   zfdum        
     76      REAL(wp) ::   zdum        
    7777      REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment 
    7878      REAL(wp) ::   zs_snic      ! snow-ice salinity 
    … …  
    9595      REAL(wp), POINTER, DIMENSION(:) ::   zq_rema     ! remaining heat at the end of the routine    (J.m-2) 
    9696      REAL(wp), POINTER, DIMENSION(:) ::   zf_tt       ! Heat budget to determine melting or freezing(W.m-2) 
     97      REAL(wp), POINTER, DIMENSION(:) ::   zevap_rema  ! remaining mass flux from sublimation        (kg.m-2) 
    9798 
    9899      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_mel   ! snow melt  
    … …  
    105106 
    106107      REAL(wp), POINTER, DIMENSION(:) ::   zqh_i       ! total ice heat content  (J.m-2) 
    107       REAL(wp), POINTER, DIMENSION(:) ::   zqh_s       ! total snow heat content (J.m-2) 
    108       REAL(wp), POINTER, DIMENSION(:) ::   zq_s        ! total snow enthalpy     (J.m-3) 
    109108      REAL(wp), POINTER, DIMENSION(:) ::   zsnw        ! distribution of snow after wind blowing 
    110109 
    … …  
    122121      END SELECT 
    123122 
    124       CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw ) 
    125       CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 
     123      CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw, zevap_rema ) 
     124      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i ) 
    126125      CALL wrk_alloc( jpij, nlay_i, zdeltah, zh_i ) 
    127126      CALL wrk_alloc( jpij, nlay_i, icount ) 
    128127        
    129       dh_i_surf  (:) = 0._wp ; dh_i_bott  (:) = 0._wp ; dh_snowice(:) = 0._wp 
     128      dh_i_surf  (:) = 0._wp ; dh_i_bott  (:) = 0._wp ; dh_snowice(:) = 0._wp ; dh_i_sub(:) = 0._wp 
    130129      dsm_i_se_1d(:) = 0._wp ; dsm_i_si_1d(:) = 0._wp    
    131130 
    132131      zqprec   (:) = 0._wp ; zq_su    (:) = 0._wp ; zq_bo    (:) = 0._wp ; zf_tt(:) = 0._wp 
    133       zq_rema  (:) = 0._wp ; zsnw     (:) = 0._wp 
     132      zq_rema  (:) = 0._wp ; zsnw     (:) = 0._wp ; zevap_rema(:) = 0._wp ; 
    134133      zdh_s_mel(:) = 0._wp ; zdh_s_pre(:) = 0._wp ; zdh_s_sub(:) = 0._wp ; zqh_i(:) = 0._wp 
    135       zqh_s    (:) = 0._wp ; zq_s     (:) = 0._wp      
    136134 
    137135      zdeltah(:,:) = 0._wp ; zh_i(:,:) = 0._wp        
    … …  
    159157      ! 
    160158      DO ji = kideb, kiut 
    161          zfdum      = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji)  
     159         zdum       = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji)  
    162160         zf_tt(ji)  = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji)  
    163161 
    164          zq_su (ji) = MAX( 0._wp, zfdum     * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) ) 
     162         zq_su (ji) = MAX( 0._wp, zdum      * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) ) 
    165163         zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice ) 
    166164      END DO 
    … …  
    187185      !  2) Computing layer thicknesses and enthalpies.            ! 
    188186      !------------------------------------------------------------! 
    189       ! 
    190       DO jk = 1, nlay_s 
    191          DO ji = kideb, kiut 
    192             zqh_s(ji) =  zqh_s(ji) + q_s_1d(ji,jk) * ht_s_1d(ji) * r1_nlay_s 
    193          END DO 
    194       END DO 
    195187      ! 
    196188      DO jk = 1, nlay_i 
    … …  
    275267      END DO 
    276268 
    277       !---------------------- 
    278       ! 3.2 Snow sublimation  
    279       !---------------------- 
     269      !------------------------------ 
     270      ! 3.2 Sublimation (part1: snow)  
     271      !------------------------------ 
    280272      ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates 
    281273      ! clem comment: not counted in mass/heat exchange in limsbc since this is an exchange with atm. (not ocean) 
    282       ! clem comment: ice should also sublimate 
    283274      zdeltah(:,:) = 0._wp 
    284       ! coupled mode: sublimation is set to 0 (evap_ice = 0) until further notice 
    285       ! forced  mode: snow thickness change due to sublimation 
    286       DO ji = kideb, kiut 
    287          zdh_s_sub(ji)  =  MAX( - ht_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice ) 
    288          ! Heat flux by sublimation [W.m-2], < 0 
    289          !      sublimate first snow that had fallen, then pre-existing snow 
     275      DO ji = kideb, kiut 
     276         zdh_s_sub(ji)  = MAX( - ht_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice ) 
     277         ! remaining evap in kg.m-2 (used for ice melting later on) 
     278         zevap_rema(ji)  = evap_ice_1d(ji) * rdt_ice + zdh_s_sub(ji) * rhosn 
     279         ! Heat flux by sublimation [W.m-2], < 0 (sublimate first snow that had fallen, then pre-existing snow) 
    290280         zdeltah(ji,1)  = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 
    291281         hfx_sub_1d(ji) = hfx_sub_1d(ji) + ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * q_s_1d(ji,1)  & 
    … …  
    309299      !------------------------------------------- 
    310300      ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow) 
    311       zq_s(:) = 0._wp  
    312301      DO jk = 1, nlay_s 
    313302         DO ji = kideb,kiut 
    314             rswitch       = MAX(  0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 )  ) 
    315             q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) *                          & 
    316               &            ( (   zdh_s_pre(ji)             ) * zqprec(ji) +  & 
    317               &              (   ht_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) ) 
    318             zq_s(ji)     =  zq_s(ji) + q_s_1d(ji,jk) 
     303            rswitch       = MAX( 0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 ) ) 
     304            q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) *           & 
     305              &            ( ( zdh_s_pre(ji)               ) * zqprec(ji) +  & 
     306              &              ( ht_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) ) 
    319307         END DO 
    320308      END DO 
    … …  
    370358               zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0] 
    371359                
    372                ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 
     360               ! Contribution to salt flux >0 (clem: using sm_i_1d and not s_i_1d(jk) is ok) 
    373361               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 
    374362                
    … …  
    383371                
    384372            END IF 
     373            ! ---------------------- 
     374            ! Sublimation part2: ice 
     375            ! ---------------------- 
     376            zdum      = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoic ) 
     377            zdeltah(ji,jk) = zdeltah(ji,jk) + zdum 
     378            dh_i_sub(ji)  = dh_i_sub(ji) + zdum 
     379            ! Salt flux > 0 (clem2016: flux is sent to the ocean for simplicity but salt should remain in the ice except if all ice is melted. 
     380            !                          It must be corrected at some point) 
     381            sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * sm_i_1d(ji) * r1_rdtice 
     382            ! Heat flux [W.m-2], < 0 
     383            hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * q_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice 
     384            ! Mass flux > 0 
     385            wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_rdtice 
     386            ! update remaining mass flux 
     387            zevap_rema(ji)  = zevap_rema(ji) + zdum * rhoic 
     388             
    385389            ! record which layers have disappeared (for bottom melting)  
    386390            !    => icount=0 : no layer has vanished 
    … …  
    389393            icount(ji,jk) = NINT( rswitch ) 
    390394            zh_i(ji,jk)   = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) ) 
    391  
     395                         
    392396            ! update heat content (J.m-2) and layer thickness 
    393397            qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk) 
    … …  
    397401      ! update ice thickness 
    398402      DO ji = kideb, kiut 
    399          ht_i_1d(ji) =  MAX( 0._wp , ht_i_1d(ji) + dh_i_surf(ji) ) 
     403         ht_i_1d(ji) =  MAX( 0._wp , ht_i_1d(ji) + dh_i_surf(ji) + dh_i_sub(ji) ) 
     404      END DO 
     405 
     406      ! remaining "potential" evap is sent to ocean 
     407      DO ji = kideb, kiut 
     408         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 
     409         wfx_err_sub(ii,ij) = wfx_err_sub(ii,ij) - zevap_rema(ji) * a_i_1d(ji) * r1_rdtice  ! <=0 (net evap for the ocean in kg.m-2.s-1) 
    400410      END DO 
    401411 
    … …  
    686696      WHERE( ht_i_1d == 0._wp ) a_i_1d = 0._wp 
    687697       
    688       CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw ) 
    689       CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 
     698      CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw, zevap_rema ) 
     699      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i ) 
    690700      CALL wrk_dealloc( jpij, nlay_i, zdeltah, zh_i ) 
    691701      CALL wrk_dealloc( jpij, nlay_i, icount ) 

  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90

    r5781 r6436  
    163163               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )   !0 if no ice and 1 if yes 
    164164               ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 
     165            END DO 
     166         END DO 
     167      END DO 
     168      ! Force the upper limit of ht_i to always be < hi_max (99 m). 
     169      DO jj = 1, jpj 
     170         DO ji = 1, jpi 
     171            rswitch = MAX( 0._wp , SIGN( 1._wp, ht_i(ji,jj,jpl) - epsi20 ) ) 
     172            ht_i(ji,jj,jpl) = MIN( ht_i(ji,jj,jpl) , hi_max(jpl) ) 
     173            a_i (ji,jj,jpl) = v_i(ji,jj,jpl) / MAX( ht_i(ji,jj,jpl) , epsi20 ) * rswitch 
     174         END DO 
     175      END DO 
     176 
     177      DO jl = 1, jpl 
     178         DO jj = 1, jpj 
     179            DO ji = 1, jpi 
     180               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )   !0 if no ice and 1 if yes 
    165181               ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 
    166182               o_i(ji,jj,jl)  = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 
    … …  
    168184         END DO 
    169185      END DO 
    170  
     186       
    171187      IF(  nn_icesal == 2  )THEN 
    172188         DO jl = 1, jpl 

  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90

    r6333 r6436  
    182182      CALL iom_put( "destrp"      , diag_trp_es         )        ! advected snw enthalpy (W/m2) 
    183183 
    184       CALL iom_put( "sfxbog"      , sfx_bog * rday      )        ! salt flux from brines 
    185       CALL iom_put( "sfxbom"      , sfx_bom * rday      )        ! salt flux from brines 
    186       CALL iom_put( "sfxsum"      , sfx_sum * rday      )        ! salt flux from brines 
    187       CALL iom_put( "sfxsni"      , sfx_sni * rday      )        ! salt flux from brines 
    188       CALL iom_put( "sfxopw"      , sfx_opw * rday      )        ! salt flux from brines 
     184      CALL iom_put( "sfxbog"      , sfx_bog * rday      )        ! salt flux from bottom growth 
     185      CALL iom_put( "sfxbom"      , sfx_bom * rday      )        ! salt flux from bottom melt 
     186      CALL iom_put( "sfxsum"      , sfx_sum * rday      )        ! salt flux from surface melt 
     187      CALL iom_put( "sfxsni"      , sfx_sni * rday      )        ! salt flux from snow ice formation 
     188      CALL iom_put( "sfxopw"      , sfx_opw * rday      )        ! salt flux from open water formation 
    189189      CALL iom_put( "sfxdyn"      , sfx_dyn * rday      )        ! salt flux from ridging rafting 
    190       CALL iom_put( "sfxres"      , sfx_res * rday      )        ! salt flux from limupdate (resultant) 
     190      CALL iom_put( "sfxres"      , sfx_res * rday      )        ! salt flux from residual 
    191191      CALL iom_put( "sfxbri"      , sfx_bri * rday      )        ! salt flux from brines 
     192      CALL iom_put( "sfxsub"      , sfx_sub * rday      )        ! salt flux from sublimation 
    192193      CALL iom_put( "sfx"         , sfx     * rday      )        ! total salt flux 
    193194 

  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90

    r6333 r6436  
    4545   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qns_ice_1d   
    4646   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   t_bo_1d      
     47   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rn_amax_1d 
    4748 
    4849   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_sum_1d 
    … …  
    5152   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_dif_1d 
    5253   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_opw_1d 
    53    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rn_amax_1d 
    5454   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_snw_1d 
    5555   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_err_1d 
    … …  
    8484   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_res_1d   
    8585 
     86   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_sub_1d 
     87 
    8688   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sprecip_1d    !: <==> the 2D  sprecip 
    8789   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   frld_1d       !: <==> the 2D  frld 
    … …  
    9294   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   evap_ice_1d   !: <==> the 2D  evap_ice 
    9395   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qprec_ice_1d  !: <==> the 2D  qprec_ice 
     96   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qevap_ice_1d  !: <==> the 3D  qevap_ice 
    9497   !                                                     ! to reintegrate longwave flux inside the ice thermodynamics 
    9598   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   i0            !: fraction of radiation transmitted to the ice 
    … …  
    108111   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_s_tot      !: Snow accretion/ablation        [m] 
    109112   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_surf     !: Ice surface accretion/ablation [m] 
     113   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_sub      !: Ice surface sublimation [m] 
    110114   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_bott     !: Ice bottom accretion/ablation  [m] 
    111115   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_snowice    !: Snow ice formation             [m of ice] 
    … …  
    155159         &      wfx_sum_1d(jpij)  , wfx_sni_1d (jpij) , wfx_opw_1d (jpij) , wfx_res_1d(jpij) ,  & 
    156160         &      dqns_ice_1d(jpij) , evap_ice_1d (jpij),                                         & 
    157          &      qprec_ice_1d(jpij), i0         (jpij) ,                                         &   
     161         &      qprec_ice_1d(jpij), qevap_ice_1d(jpij), i0         (jpij) ,                     &   
    158162         &      sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) , sfx_sum_1d (jpij),  & 
    159          &      sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) ,                     & 
     163         &      sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , sfx_sub_1d (jpij),  & 
    160164         &      dsm_i_fl_1d(jpij) , dsm_i_gd_1d(jpij) , dsm_i_se_1d(jpij) ,                     &      
    161165         &      dsm_i_si_1d(jpij) , hicol_1d    (jpij)                     , STAT=ierr(2) ) 
    … …  
    163167      ALLOCATE( t_su_1d   (jpij) , a_i_1d   (jpij) , ht_i_1d  (jpij) ,    &    
    164168         &      ht_s_1d   (jpij) , fc_su    (jpij) , fc_bo_i  (jpij) ,    &     
    165          &      dh_s_tot  (jpij) , dh_i_surf(jpij) , dh_i_bott(jpij) ,    &     
    166          &      dh_snowice(jpij) , sm_i_1d  (jpij) , s_i_new  (jpij) ,    & 
     169         &      dh_s_tot  (jpij) , dh_i_surf(jpij) , dh_i_sub (jpij) ,    &     
     170         &      dh_i_bott (jpij) ,dh_snowice(jpij) , sm_i_1d  (jpij) , s_i_new  (jpij) ,    & 
    167171         &      t_s_1d(jpij,nlay_s) , t_i_1d(jpij,nlay_i) , s_i_1d(jpij,nlay_i) ,  &             
    168172         &      q_i_1d(jpij,nlay_i+1) , q_s_1d(jpij,nlay_s) ,                        & 

  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r6333 r6436  
    150150      CALL iom_put("e3v_0", e3t_0(:,:,:) ) 
    151151      ! 
    152       IF( .NOT.lk_vvl ) THEN 
    153          CALL iom_put( "e3t" , fse3t_n(:,:,:) ) 
    154          CALL iom_put( "e3u" , fse3u_n(:,:,:) ) 
    155          CALL iom_put( "e3v" , fse3v_n(:,:,:) ) 
    156          CALL iom_put( "e3w" , fse3w_n(:,:,:) ) 
    157       ENDIF 
     152      CALL iom_put( "e3t" , fse3t_n(:,:,:) ) 
     153      CALL iom_put( "e3u" , fse3u_n(:,:,:) ) 
     154      CALL iom_put( "e3v" , fse3v_n(:,:,:) ) 
     155      CALL iom_put( "e3w" , fse3w_n(:,:,:) ) 
     156      IF( iom_use("e3tdef") )   & 
     157         CALL iom_put( "e3tdef"  , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 
     158 
    158159 
    159160      CALL iom_put( "ssh" , sshn )                 ! sea surface height 
    … …  
    247248      CALL iom_put( "avm" , avmu                       )    ! T vert. eddy visc. coef. 
    248249      CALL iom_put( "avs" , fsavs(:,:,:)               )    ! S vert. eddy diff. coef. (useful only with key_zdfddm) 
     250                                                            ! Log of eddy diff coef 
     251      IF( iom_use('logavt') )   CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt  (:,:,:) ) ) ) 
     252      IF( iom_use('logavs') )   CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, fsavs(:,:,:) ) ) ) 
    249253 
    250254      IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN 
    … …  
    311315         CALL iom_put( "eken", rke )            
    312316      ENDIF 
    313           
     317      ! 
     318      CALL iom_put( "hdiv", hdivn )                  ! Horizontal divergence 
     319      ! 
    314320      IF( iom_use("u_masstr") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN 
    315321         z3d(:,:,jpk) = 0.e0 

  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r5781 r6436  
    665665         ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 
    666666      END DO 
    667  
    668       ! Write outputs 
    669       ! ============= 
    670       CALL iom_put(     "e3t" , fse3t_n  (:,:,:) ) 
    671       CALL iom_put(     "e3u" , fse3u_n  (:,:,:) ) 
    672       CALL iom_put(     "e3v" , fse3v_n  (:,:,:) ) 
    673       CALL iom_put(     "e3w" , fse3w_n  (:,:,:) ) 
    674       CALL iom_put( "tpt_dep" , fsde3w_n (:,:,:) ) 
    675       IF( iom_use("e3tdef") )   & 
    676          CALL iom_put( "e3tdef"  , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 
    677667 
    678668      ! write restart file 

  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/LBC/mppini.F90

    r5783 r6436  
    201201       
    202202#endif 
    203       IF(lwp) THEN 
    204          WRITE(numout,*) 
    205          WRITE(numout,*) '           defines mpp subdomains' 
    206          WRITE(numout,*) '           ----------------------' 
    207          WRITE(numout,*) '           iresti=',iresti,' irestj=',irestj 
    208          WRITE(numout,*) '           jpni  =',jpni  ,' jpnj  =',jpnj 
    209          ifreq = 4 
    210          il1   = 1 
    211          DO jn = 1, (jpni-1)/ifreq+1 
    212             il2 = MIN( jpni, il1+ifreq-1 ) 
    213             WRITE(numout,*) 
    214             WRITE(numout,9200) ('***',ji = il1,il2-1) 
    215             DO jj = jpnj, 1, -1 
    216                WRITE(numout,9203) ('   ',ji = il1,il2-1) 
    217                WRITE(numout,9202) jj, ( ilcit(ji,jj),ilcjt(ji,jj),ji = il1,il2 ) 
    218                WRITE(numout,9203) ('   ',ji = il1,il2-1) 
    219                WRITE(numout,9200) ('***',ji = il1,il2-1) 
    220             END DO 
    221             WRITE(numout,9201) (ji,ji = il1,il2) 
    222             il1 = il1+ifreq 
    223          END DO 
    224  9200    FORMAT('     ***',20('*************',a3)) 
    225  9203    FORMAT('     *     ',20('         *   ',a3)) 
    226  9201    FORMAT('        ',20('   ',i3,'          ')) 
    227  9202    FORMAT(' ',i3,' *  ',20(i3,'  x',i3,'   *   ')) 
    228       ENDIF 
    229  
    230       zidom = nreci 
    231       DO ji = 1, jpni 
    232          zidom = zidom + ilcit(ji,1) - nreci 
    233       END DO 
    234       IF(lwp) WRITE(numout,*) 
    235       IF(lwp) WRITE(numout,*)' sum ilcit(i,1) = ', zidom, ' jpiglo = ', jpiglo 
    236        
    237       zjdom = nrecj 
    238       DO jj = 1, jpnj 
    239          zjdom = zjdom + ilcjt(1,jj) - nrecj 
    240       END DO 
    241       IF(lwp) WRITE(numout,*)' sum ilcit(1,j) = ', zjdom, ' jpjglo = ', jpjglo 
    242       IF(lwp) WRITE(numout,*) 
    243        
    244203 
    245204      !  2. Index arrays for subdomains 
    … …  
    304263         nlejt(jn) = nlej 
    305264      END DO 
    306        
    307  
    308       ! 4. From global to local 
     265 
     266      ! 4. Subdomain print 
     267      ! ------------------ 
     268       
     269      IF(lwp) WRITE(numout,*) 
     270      IF(lwp) WRITE(numout,*) ' mpp_init: defines mpp subdomains' 
     271      IF(lwp) WRITE(numout,*) ' ~~~~~~  ----------------------' 
     272      IF(lwp) WRITE(numout,*) 
     273      IF(lwp) WRITE(numout,*) 'iresti=',iresti,' irestj=',irestj 
     274      IF(lwp) WRITE(numout,*) 
     275      IF(lwp) WRITE(numout,*) 'jpni=',jpni,' jpnj=',jpnj 
     276      zidom = nreci 
     277      DO ji = 1, jpni 
     278         zidom = zidom + ilcit(ji,1) - nreci 
     279      END DO 
     280      IF(lwp) WRITE(numout,*) 
     281      IF(lwp) WRITE(numout,*)' sum ilcit(i,1)=', zidom, ' jpiglo=', jpiglo 
     282 
     283      zjdom = nrecj 
     284      DO jj = 1, jpnj 
     285         zjdom = zjdom + ilcjt(1,jj) - nrecj 
     286      END DO 
     287      IF(lwp) WRITE(numout,*)' sum ilcit(1,j)=', zjdom, ' jpjglo=', jpjglo 
     288      IF(lwp) WRITE(numout,*) 
     289 
     290      IF(lwp) THEN 
     291         ifreq = 4 
     292         il1   = 1 
     293         DO jn = 1, (jpni-1)/ifreq+1 
     294            il2 = MIN( jpni, il1+ifreq-1 ) 
     295            WRITE(numout,*) 
     296            WRITE(numout,9200) ('***',ji = il1,il2-1) 
     297            DO jj = jpnj, 1, -1 
     298               WRITE(numout,9203) ('   ',ji = il1,il2-1) 
     299               WRITE(numout,9202) jj, ( ilcit(ji,jj),ilcjt(ji,jj),ji = il1,il2 ) 
     300               WRITE(numout,9204) (nfipproc(ji,jj),ji=il1,il2) 
     301               WRITE(numout,9203) ('   ',ji = il1,il2-1) 
     302               WRITE(numout,9200) ('***',ji = il1,il2-1) 
     303            END DO 
     304            WRITE(numout,9201) (ji,ji = il1,il2) 
     305            il1 = il1+ifreq 
     306         END DO 
     307 9200     FORMAT('     ***',20('*************',a3)) 
     308 9203     FORMAT('     *     ',20('         *   ',a3)) 
     309 9201     FORMAT('        ',20('   ',i3,'          ')) 
     310 9202     FORMAT(' ',i3,' *  ',20(i3,'  x',i3,'   *   ')) 
     311 9204     FORMAT('     *  ',20('      ',i3,'   *   ')) 
     312      ENDIF 
     313 
     314      ! 5. From global to local 
    309315      ! ----------------------- 
    310316 
    … …  
    313319 
    314320 
    315       ! 5. Subdomain neighbours 
     321      ! 6. Subdomain neighbours 
    316322      ! ---------------------- 
    317323 
    … …  
    436442         WRITE(numout,*) ' nimpp  = ', nimpp 
    437443         WRITE(numout,*) ' njmpp  = ', njmpp 
    438          WRITE(numout,*) ' nbse   = ', nbse  , ' npse   = ', npse 
    439          WRITE(numout,*) ' nbsw   = ', nbsw  , ' npsw   = ', npsw 
    440          WRITE(numout,*) ' nbne   = ', nbne  , ' npne   = ', npne 
    441          WRITE(numout,*) ' nbnw   = ', nbnw  , ' npnw   = ', npnw 
     444         WRITE(numout,*) ' nreci  = ', nreci  , ' npse   = ', npse 
     445         WRITE(numout,*) ' nrecj  = ', nrecj  , ' npsw   = ', npsw 
     446         WRITE(numout,*) ' jpreci = ', jpreci , ' npne   = ', npne 
     447         WRITE(numout,*) ' jprecj = ', jprecj , ' npnw   = ', npnw 
     448         WRITE(numout,*) 
    442449      ENDIF 
    443450 
    … …  
    446453      ! Prepare mpp north fold 
    447454 
    448       IF (jperio >= 3 .AND. jperio <= 6 .AND. jpni > 1 ) THEN 
     455      IF( jperio >= 3 .AND. jperio <= 6 .AND. jpni > 1 ) THEN 
    449456         CALL mpp_ini_north 
    450       END IF 
     457         IF(lwp) WRITE(numout,*) ' mpp_init : North fold boundary prepared for jpni >1' 
     458      ENDIF 
    451459 
    452460      ! Prepare NetCDF output file (if necessary) 

  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/LBC/mppini_2.h90

    r5781 r6436  
    318318         ENDIF 
    319319 
     320         ! Check wet points over the entire domain to preserve the MPI communication stencil 
    320321         isurf = 0 
    321          DO jj = 1+jprecj, ilj-jprecj 
    322             DO  ji = 1+jpreci, ili-jpreci 
     322         DO jj = 1, ilj 
     323            DO  ji = 1, ili 
    323324               IF( imask(ji+iimppt(ii,ij)-1, jj+ijmppt(ii,ij)-1) == 1) isurf = isurf+1 
    324325            END DO 
    325326         END DO 
     327 
    326328         IF(isurf /= 0) THEN 
    327329            icont = icont + 1 
    … …  
    333335 
    334336      nfipproc(:,:) = ipproc(:,:) 
    335  
    336337 
    337338      ! Control 
    … …  
    441442      ii = iin(narea) 
    442443      ij = ijn(narea) 
     444 
     445      ! set default neighbours 
     446      noso = ioso(ii,ij) 
     447      nowe = iowe(ii,ij) 
     448      noea = ioea(ii,ij) 
     449      nono = iono(ii,ij)  
     450      npse = iose(ii,ij) 
     451      npsw = iosw(ii,ij) 
     452      npne = ione(ii,ij) 
     453      npnw = ionw(ii,ij) 
     454 
     455      ! check neighbours location 
    443456      IF( ioso(ii,ij) >= 0 .AND. ioso(ii,ij) <= (jpni*jpnj-1) ) THEN  
    444457         iiso = 1 + MOD(ioso(ii,ij),jpni) 
    … …  
    511524      IF (lwp) THEN 
    512525         CALL ctl_opn( inum, 'layout.dat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE., narea ) 
     526         WRITE(inum,'(a)') '   jpnij     jpi     jpj     jpk  jpiglo  jpjglo' 
    513527         WRITE(inum,'(6i8)') jpnij,jpi,jpj,jpk,jpiglo,jpjglo 
    514528         WRITE(inum,'(a)') 'NAREA nlci nlcj nldi nldj nlei nlej nimpp njmpp' 
    … …  
    523537      END IF 
    524538 
    525       IF( nperio == 1 .AND.jpni /= 1 ) CALL ctl_stop( ' mpp_init2:  error on cyclicity' ) 
    526  
    527       ! Prepare mpp north fold 
    528  
    529       IF( jperio >= 3 .AND. jperio <= 6 .AND. jpni > 1 ) THEN 
    530          CALL mpp_ini_north 
    531          IF(lwp) WRITE(numout,*) ' mpp_init2 : North fold boundary prepared for jpni >1' 
    532       ENDIF 
    533  
    534539      ! Defined npolj, either 0, 3 , 4 , 5 , 6 
    535540      ! In this case the important thing is that npolj /= 0 
    … …  
    548553      ENDIF 
    549554 
     555      ! Periodicity : no corner if nbondi = 2 and nperio != 1 
     556 
     557      IF(lwp) THEN 
     558         WRITE(numout,*) ' nproc  = ', nproc 
     559         WRITE(numout,*) ' nowe   = ', nowe  , ' noea   =  ', noea 
     560         WRITE(numout,*) ' nono   = ', nono  , ' noso   =  ', noso 
     561         WRITE(numout,*) ' nbondi = ', nbondi 
     562         WRITE(numout,*) ' nbondj = ', nbondj 
     563         WRITE(numout,*) ' npolj  = ', npolj 
     564         WRITE(numout,*) ' nperio = ', nperio 
     565         WRITE(numout,*) ' nlci   = ', nlci 
     566         WRITE(numout,*) ' nlcj   = ', nlcj 
     567         WRITE(numout,*) ' nimpp  = ', nimpp 
     568         WRITE(numout,*) ' njmpp  = ', njmpp 
     569         WRITE(numout,*) ' nreci  = ', nreci  , ' npse   = ', npse 
     570         WRITE(numout,*) ' nrecj  = ', nrecj  , ' npsw   = ', npsw 
     571         WRITE(numout,*) ' jpreci = ', jpreci , ' npne   = ', npne 
     572         WRITE(numout,*) ' jprecj = ', jprecj , ' npnw   = ', npnw 
     573         WRITE(numout,*) 
     574      ENDIF 
     575 
     576      IF( nperio == 1 .AND. jpni /= 1 ) CALL ctl_stop( ' mpp_init2: error on cyclicity' ) 
     577 
     578      ! Prepare mpp north fold 
     579 
     580      IF( jperio >= 3 .AND. jperio <= 6 .AND. jpni > 1 ) THEN 
     581         CALL mpp_ini_north 
     582         IF(lwp) WRITE(numout,*) ' mpp_init2 : North fold boundary prepared for jpni >1' 
     583      ENDIF 
     584 
    550585      ! Prepare NetCDF output file (if necessary) 
    551586      CALL mpp_init_ioipsl 
    552587 
    553       ! Periodicity : no corner if nbondi = 2 and nperio != 1 
    554  
    555       IF(lwp) THEN 
    556          WRITE(numout,*) ' nproc=  ',nproc 
    557          WRITE(numout,*) ' nowe=   ',nowe 
    558          WRITE(numout,*) ' noea=   ',noea 
    559          WRITE(numout,*) ' nono=   ',nono 
    560          WRITE(numout,*) ' noso=   ',noso 
    561          WRITE(numout,*) ' nbondi= ',nbondi 
    562          WRITE(numout,*) ' nbondj= ',nbondj 
    563          WRITE(numout,*) ' npolj=  ',npolj 
    564          WRITE(numout,*) ' nperio= ',nperio 
    565          WRITE(numout,*) ' nlci=   ',nlci 
    566          WRITE(numout,*) ' nlcj=   ',nlcj 
    567          WRITE(numout,*) ' nimpp=  ',nimpp 
    568          WRITE(numout,*) ' njmpp=  ',njmpp 
    569          WRITE(numout,*) ' nbse=   ',nbse,' npse= ',npse 
    570          WRITE(numout,*) ' nbsw=   ',nbsw,' npsw= ',npsw 
    571          WRITE(numout,*) ' nbne=   ',nbne,' npne= ',npne 
    572          WRITE(numout,*) ' nbnw=   ',nbnw,' npnw= ',npnw 
    573       ENDIF 
    574588 
    575589   END SUBROUTINE mpp_init2 

  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90

    r5781 r6436  
    188188            DO jj = 2, jpjm1 
    189189               DO ji = fs_2, fs_jpim1   ! vector opt. 
    190                   IF (miku(ji,jj) .GT. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji  ,jj  ),                   5._wp) 
    191                   IF (miku(ji,jj) .LT. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji+1,jj  ),                   5._wp) 
    192                   IF (miku(ji,jj) .EQ. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji  ,jj  ), hmlpt(ji+1,jj  ), 5._wp) 
    193                   IF (mikv(ji,jj) .GT. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj  ),                   5._wp) 
    194                   IF (mikv(ji,jj) .LT. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj+1),                   5._wp) 
    195                   IF (mikv(ji,jj) .EQ. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj  ), hmlpt(ji  ,jj+1), 5._wp) 
     190               zhmlpu(ji,jj) = ( MAX(hmlpt(ji,jj)  , hmlpt  (ji+1,jj  ), 5._wp)   & 
     191                  &            - MAX(risfdep(ji,jj), risfdep(ji+1,jj  )       )   ) 
     192               zhmlpv(ji,jj) = ( MAX(hmlpt  (ji,jj), hmlpt  (ji  ,jj+1), 5._wp)   & 
     193                  &            - MAX(risfdep(ji,jj), risfdep(ji  ,jj+1)       )   ) 
    196194               ENDDO 
    197195            ENDDO 

  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90

    r5783 r6436  
    8080   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qemp_oce       !: heat flux of precip and evap over ocean     [W/m2] 
    8181   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qemp_ice       !: heat flux of precip and evap over ice       [W/m2] 
    82    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qprec_ice      !: heat flux of precip over ice                [J/m3] 
     82   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qevap_ice      !: heat flux of evap over ice                  [W/m2] 
     83   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qprec_ice      !: enthalpy of precip over ice                 [J/m3] 
    8384   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   emp_oce        !: evap - precip over ocean                 [kg/m2/s] 
    8485#endif 
    … …  
    144145#endif 
    145146#if defined key_lim3 
    146          &      evap_ice(jpi,jpj,jpl) , devap_ice(jpi,jpj,jpl) , qprec_ice(jpi,jpj) ,  & 
    147          &      qemp_ice(jpi,jpj)     , qemp_oce(jpi,jpj)      ,                       & 
    148          &      qns_oce (jpi,jpj)     , qsr_oce (jpi,jpj)      , emp_oce (jpi,jpj)  ,  & 
     147         &      evap_ice(jpi,jpj,jpl) , devap_ice(jpi,jpj,jpl) , qprec_ice(jpi,jpj) ,   & 
     148         &      qemp_ice(jpi,jpj)     , qevap_ice(jpi,jpj,jpl) , qemp_oce (jpi,jpj) ,   & 
     149         &      qns_oce (jpi,jpj)     , qsr_oce  (jpi,jpj)     , emp_oce (jpi,jpj)  ,   & 
    149150#endif 
    150151         &      emp_ice(jpi,jpj)      ,  STAT= ierr(1) ) 

  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90

    r5783 r6436  
    684684      qprec_ice(:,:) = rhosn * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 
    685685 
     686      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- ! 
     687      DO jl = 1, jpl 
     688         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) - lfus ) 
     689                                   ! but then qemp_ice should also include sublimation  
     690      END DO 
     691 
    686692      CALL wrk_dealloc( jpi,jpj, zevap, zsnw )  
    687693#endif 

  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r5783 r6436  
    612612      ! --- evaporation --- ! 
    613613      z1_lsub = 1._wp / Lsub 
    614       evap_ice (:,:,:) = qla_ice (:,:,:) * z1_lsub ! sublimation 
    615       devap_ice(:,:,:) = dqla_ice(:,:,:) * z1_lsub 
    616       zevap    (:,:)   = emp(:,:) + tprecip(:,:)   ! evaporation over ocean 
     614      evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_lsub    ! sublimation 
     615      devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_lsub    ! d(sublimation)/dT 
     616      zevap    (:,:)   = rn_efac * ( emp(:,:) + tprecip(:,:) )  ! evaporation over ocean 
    617617 
    618618      ! --- evaporation minus precipitation --- ! 
    … …  
    637637      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- ! 
    638638      qprec_ice(:,:) = rhosn * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 
     639 
     640      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- ! 
     641      DO jl = 1, jpl 
     642         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) ) 
     643                                   ! But we do not have Tice => consider it at 0°C => evap=0  
     644      END DO 
    639645 
    640646      CALL wrk_dealloc( jpi,jpj, zevap, zsnw )  

  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r6237 r6436  
    13781378      ! 
    13791379      INTEGER ::   jl         ! dummy loop index 
    1380       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zcptn, ztmp, zicefr, zmsk 
    1381       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot 
    1382       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice 
    1383       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zevap, zsnw, zqns_oce, zqsr_oce, zqprec_ice, zqemp_oce ! for LIM3 
     1380      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 
     1382      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice 
     1383      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice 
    13841384      !!---------------------------------------------------------------------- 
    13851385      ! 
    13861386      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_flx') 
    13871387      ! 
    1388       CALL wrk_alloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot ) 
    1389       CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice ) 
     1388      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 ) 
     1390      CALL wrk_alloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice ) 
     1391      CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice ) 
    13901392 
    13911393      IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0) 
    … …  
    14231425      END SELECT 
    14241426 
    1425       IF( iom_use('subl_ai_cea') )   & 
    1426          CALL iom_put( 'subl_ai_cea', frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) )   ! Sublimation over sea-ice         (cell average) 
    1427       !    
    1428       !                                                           ! runoffs and calving (put in emp_tot) 
     1427#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 ) 
     1431       
     1432      ! --- evaporation (kg/m2/s) --- ! 
     1433      zevap_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1) 
     1434      ! since the sensitivity of evap to temperature (devap/dT) is not prescribed by the atmosphere, we set it to 0 
     1435      ! therefore, sublimation is not redistributed over the ice categories in case no subgrid scale fluxes are provided by atm. 
     1436      zdevap_ice(:,:) = 0._wp 
     1437       
     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) 
     1445      IF( srcv(jpr_rnf)%laction )   rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1) 
     1446      IF( srcv(jpr_cal)%laction ) THEN  
     1447         zemp_tot(:,:) = zemp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1) 
     1448         CALL iom_put( 'calving_cea', frcv(jpr_cal)%z3(:,:,1) ) 
     1449      ENDIF 
     1450 
     1451      IF( ln_mixcpl ) THEN 
     1452         emp_tot(:,:) = emp_tot(:,:) * xcplmask(:,:,0) + zemp_tot(:,:) * zmsk(:,:) 
     1453         emp_ice(:,:) = emp_ice(:,:) * xcplmask(:,:,0) + zemp_ice(:,:) * zmsk(:,:) 
     1454         emp_oce(:,:) = emp_oce(:,:) * xcplmask(:,:,0) + zemp_oce(:,:) * zmsk(:,:) 
     1455         sprecip(:,:) = sprecip(:,:) * xcplmask(:,:,0) + zsprecip(:,:) * zmsk(:,:) 
     1456         tprecip(:,:) = tprecip(:,:) * xcplmask(:,:,0) + ztprecip(:,:) * zmsk(:,:) 
     1457         DO jl=1,jpl 
     1458            evap_ice (:,:,jl) = evap_ice (:,:,jl) * xcplmask(:,:,0) + zevap_ice (:,:) * zmsk(:,:) 
     1459            devap_ice(:,:,jl) = devap_ice(:,:,jl) * xcplmask(:,:,0) + zdevap_ice(:,:) * zmsk(:,:) 
     1460         ENDDO 
     1461      ELSE 
     1462         emp_tot(:,:) =         zemp_tot(:,:) 
     1463         emp_ice(:,:) =         zemp_ice(:,:) 
     1464         emp_oce(:,:) =         zemp_oce(:,:)      
     1465         sprecip(:,:) =         zsprecip(:,:) 
     1466         tprecip(:,:) =         ztprecip(:,:) 
     1467         DO jl=1,jpl 
     1468            evap_ice (:,:,jl) = zevap_ice (:,:) 
     1469            devap_ice(:,:,jl) = zdevap_ice(:,:) 
     1470         ENDDO 
     1471      ENDIF 
     1472 
     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)     
     1476#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(:,:) ) 
     1479      ! runoffs and calving (put in emp_tot) 
    14291480      IF( srcv(jpr_rnf)%laction )   rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1) 
    14301481      IF( srcv(jpr_cal)%laction ) THEN  
    … …  
    14501501      IF( iom_use('snow_ai_cea') )   & 
    14511502         CALL iom_put( 'snow_ai_cea', sprecip(:,:) * zicefr(:,:)             )   ! Snow        over sea-ice         (cell average) 
     1503#endif 
    14521504 
    14531505      !                                                      ! ========================= ! 
    … …  
    15051557      IF( iom_use('hflx_snow_cea') )    CALL iom_put( 'hflx_snow_cea', ztmp + sprecip(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average) 
    15061558 
    1507 #if defined key_lim3 
    1508       CALL wrk_alloc( jpi,jpj, zevap, zsnw, zqns_oce, zqprec_ice, zqemp_oce )  
    1509  
     1559#if defined key_lim3       
    15101560      ! --- evaporation --- ! 
    1511       ! clem: evap_ice is set to 0 for LIM3 since we still do not know what to do with sublimation 
    1512       ! the problem is: the atm. imposes both mass evaporation and heat removed from the snow/ice 
    1513       !                 but it is incoherent WITH the ice model   
    1514       DO jl=1,jpl 
    1515          evap_ice(:,:,jl) = 0._wp  ! should be: frcv(jpr_ievp)%z3(:,:,1) 
    1516       ENDDO 
    15171561      zevap(:,:) = zemp_tot(:,:) + ztprecip(:,:) ! evaporation over ocean 
    1518  
    1519       ! --- evaporation minus precipitation --- ! 
    1520       emp_oce(:,:) = emp_tot(:,:) - emp_ice(:,:) 
    15211562 
    15221563      ! --- non solar flux over ocean --- ! 
    … …  
    15251566      WHERE( p_frld /= 0._wp )  zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / p_frld(:,:) 
    15261567 
    1527       ! --- heat flux associated with emp --- ! 
    1528       zsnw(:,:) = 0._wp 
    1529       CALL lim_thd_snwblow( p_frld, zsnw )  ! snow distribution over ice after wind blowing 
     1568      ! --- heat flux associated with emp (W/m2) --- ! 
    15301569      zqemp_oce(:,:) = -      zevap(:,:)                   * p_frld(:,:)      *   zcptn(:,:)   &      ! evap 
    15311570         &             + ( ztprecip(:,:) - zsprecip(:,:) )                    *   zcptn(:,:)   &      ! liquid precip 
    15321571         &             +   zsprecip(:,:)                   * ( 1._wp - zsnw ) * ( zcptn(:,:) - lfus ) ! solid precip over ocean 
    1533       qemp_ice(:,:)  = -   frcv(jpr_ievp)%z3(:,:,1)        * zicefr(:,:)      *   zcptn(:,:)   &      ! ice evap 
    1534          &             +   zsprecip(:,:)                   * zsnw             * ( zcptn(:,:) - lfus ) ! solid precip over ice 
    1535  
     1572!      zqemp_ice(:,:) = -   frcv(jpr_ievp)%z3(:,:,1)        * zicefr(:,:)      *   zcptn(:,:)   &      ! ice evap 
     1573!         &             +   zsprecip(:,:)                   * zsnw             * ( zcptn(:,:) - lfus ) ! solid precip over ice 
     1574      zqemp_ice(:,:) =      zsprecip(:,:)                   * zsnw             * ( zcptn(:,:) - lfus ) ! solid precip over ice (only) 
     1575                                                                                                       ! qevap_ice=0 since we consider Tice=0°C 
     1576       
    15361577      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- ! 
    15371578      zqprec_ice(:,:) = rhosn * ( zcptn(:,:) - lfus ) 
    15381579 
    1539       ! --- total non solar flux --- ! 
    1540       zqns_tot(:,:) = zqns_tot(:,:) + qemp_ice(:,:) + zqemp_oce(:,:) 
     1580      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- ! 
     1581      DO jl = 1, jpl 
     1582         zqevap_ice(:,:,jl) = 0._wp ! should be -evap * ( ( Tice - rt0 ) * cpic ) but we do not have Tice, so we consider Tice=0°C 
     1583      END DO 
     1584 
     1585      ! --- total non solar flux (including evap/precip) --- ! 
     1586      zqns_tot(:,:) = zqns_tot(:,:) + zqemp_ice(:,:) + zqemp_oce(:,:) 
    15411587 
    15421588      ! --- in case both coupled/forced are active, we must mix values --- !  
    … …  
    15451591         qns_oce(:,:) = qns_oce(:,:) * xcplmask(:,:,0) + zqns_oce(:,:)* zmsk(:,:) 
    15461592         DO jl=1,jpl 
    1547             qns_ice(:,:,jl) = qns_ice(:,:,jl) * xcplmask(:,:,0) +  zqns_ice(:,:,jl)* zmsk(:,:) 
     1593            qns_ice  (:,:,jl) = qns_ice  (:,:,jl) * xcplmask(:,:,0) +  zqns_ice  (:,:,jl)* zmsk(:,:) 
     1594            qevap_ice(:,:,jl) = qevap_ice(:,:,jl) * xcplmask(:,:,0) +  zqevap_ice(:,:,jl)* zmsk(:,:) 
    15481595         ENDDO 
    15491596         qprec_ice(:,:) = qprec_ice(:,:) * xcplmask(:,:,0) + zqprec_ice(:,:)* zmsk(:,:) 
    15501597         qemp_oce (:,:) =  qemp_oce(:,:) * xcplmask(:,:,0) +  zqemp_oce(:,:)* zmsk(:,:) 
    1551 !!clem         evap_ice(:,:) = evap_ice(:,:) * xcplmask(:,:,0) 
     1598         qemp_ice (:,:) =  qemp_ice(:,:) * xcplmask(:,:,0) +  zqemp_ice(:,:)* zmsk(:,:) 
    15521599      ELSE 
    15531600         qns_tot  (:,:  ) = zqns_tot  (:,:  ) 
    15541601         qns_oce  (:,:  ) = zqns_oce  (:,:  ) 
    15551602         qns_ice  (:,:,:) = zqns_ice  (:,:,:) 
    1556          qprec_ice(:,:)   = zqprec_ice(:,:) 
    1557          qemp_oce (:,:)   = zqemp_oce (:,:) 
    1558       ENDIF 
    1559  
    1560       CALL wrk_dealloc( jpi,jpj, zevap, zsnw, zqns_oce, zqprec_ice, zqemp_oce )  
     1603         qevap_ice(:,:,:) = zqevap_ice(:,:,:) 
     1604         qprec_ice(:,:  ) = zqprec_ice(:,:  ) 
     1605         qemp_oce (:,:  ) = zqemp_oce (:,:  ) 
     1606         qemp_ice (:,:  ) = zqemp_ice (:,:  ) 
     1607      ENDIF 
    15611608#else 
    1562  
    15631609      ! clem: this formulation is certainly wrong... but better than it was... 
    15641610      zqns_tot(:,:) = zqns_tot(:,:)                       &            ! zqns_tot update over free ocean with: 
    … …  
    15771623         qns_ice(:,:,:) = zqns_ice(:,:,:) 
    15781624      ENDIF 
    1579  
    15801625#endif 
    15811626 
    … …  
    16281673 
    16291674#if defined key_lim3 
    1630       CALL wrk_alloc( jpi,jpj, zqsr_oce )  
    16311675      ! --- solar flux over ocean --- ! 
    16321676      !         note: p_frld cannot be = 0 since we limit the ice concentration to amax 
    … …  
    16361680      IF( ln_mixcpl ) THEN   ;   qsr_oce(:,:) = qsr_oce(:,:) * xcplmask(:,:,0) +  zqsr_oce(:,:)* zmsk(:,:) 
    16371681      ELSE                   ;   qsr_oce(:,:) = zqsr_oce(:,:)   ;   ENDIF 
    1638  
    1639       CALL wrk_dealloc( jpi,jpj, zqsr_oce )  
    16401682#endif 
    16411683 
    … …  
    16881730      fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
    16891731 
    1690       CALL wrk_dealloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot ) 
    1691       CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice ) 
     1732      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 ) 
     1734      CALL wrk_dealloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice ) 
     1735      CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice ) 
    16921736      ! 
    16931737      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_ice_flx') 

  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r6333 r6436  
    203203            ! In CLIO the cloud fraction is read in the climatology and the all-sky albedo  
    204204            ! (alb_ice) is computed within the bulk routine 
    205             CALL blk_ice_clio_flx( t_su, zalb_cs, zalb_os, alb_ice ) 
     205                                 CALL blk_ice_clio_flx( t_su, zalb_cs, zalb_os, alb_ice ) 
    206206            IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su ) 
    207207            IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
    … …  
    209209            ! albedo depends on cloud fraction because of non-linear spectral effects 
    210210            alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
    211             CALL blk_ice_core_flx( t_su, alb_ice ) 
     211                                 CALL blk_ice_core_flx( t_su, alb_ice ) 
    212212            IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su ) 
    213213            IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
    … …  
    216216            alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
    217217                                 CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su ) 
    218             ! clem: evap_ice is forced to 0 in coupled mode for now  
    219             !       but it needs to be changed (along with modif in limthd_dh) once heat flux from evap will be avail. from atm. models 
    220             evap_ice  (:,:,:) = 0._wp   ;   devap_ice (:,:,:) = 0._wp 
    221218            IF( nn_limflx == 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
    222219         END SELECT 
    … …  
    588585      sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp 
    589586      sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp 
    590       sfx_res(:,:) = 0._wp 
     587      sfx_res(:,:) = 0._wp   ;   sfx_sub(:,:) = 0._wp 
    591588       
    592589      wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp 
    … …  
    604601      hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp  
    605602      hfx_err(:,:) = 0._wp   ;   hfx_err_rem(:,:) = 0._wp 
    606       hfx_err_dif(:,:) = 0._wp   ; 
    607  
     603      hfx_err_dif(:,:) = 0._wp 
     604      wfx_err_sub(:,:) = 0._wp 
     605       
    608606      afx_tot(:,:) = 0._wp   ; 
    609607      afx_dyn(:,:) = 0._wp   ;   afx_thd(:,:) = 0._wp 

  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r5783 r6436  
    456456      !                                                ! ---------------------------------------- ! 
    457457      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 
    458          CALL iom_put( "empmr" , emp  - rnf )                   ! upward water flux 
     458         CALL iom_put( "empmr"  , emp    - rnf )                ! upward water flux 
     459         CALL iom_put( "empbmr" , emp_b  - rnf )                ! before upward water flux ( needed to recalculate the time evolution of ssh in offline ) 
    459460         CALL iom_put( "saltflx", sfx  )                        ! downward salt flux   
    460461                                                                ! (includes virtual salt flux beneath ice  

  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90

    r6333 r6436  
    215215      IF( ierr == 1 )   CALL ctl_stop( ' iso-level in z-coordinate - partial step, not allowed' ) 
    216216      IF( ierr == 2 )   CALL ctl_stop( ' isoneutral bilaplacian operator does not exist' ) 
     217      IF( ln_traldf_grif .AND. ln_isfcav         )   & 
     218           CALL ctl_stop( ' ice shelf and traldf_grif not tested') 
    217219      IF( lk_traldf_eiv .AND. .NOT.ln_traldf_iso )   & 
    218220           CALL ctl_stop( '          eddy induced velocity on tracers',   & 

  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90

    r6237 r6436  
    8080      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    8181      ! 
    82       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    83       INTEGER  ::   iikn, iiki, ikt, imkt  ! local integer 
    84       REAL(wp) ::   zN2_c        ! local scalar 
     82      INTEGER  ::   ji, jj, jk      ! dummy loop indices 
     83      INTEGER  ::   iikn, iiki, ikt ! local integer 
     84      REAL(wp) ::   zN2_c           ! local scalar 
    8585      INTEGER, POINTER, DIMENSION(:,:) ::   imld   ! 2D workspace 
    8686      !!---------------------------------------------------------------------- 
    … …  
    117117         DO jj = 1, jpj 
    118118            DO ji = 1, jpi 
    119                imkt = mikt(ji,jj) 
    120                IF( avt (ji,jj,jk) < avt_c )   imld(ji,jj) = MAX( imkt, jk )      ! Turbocline  
     119               IF( avt (ji,jj,jk) < avt_c * wmask(ji,jj,jk) )   imld(ji,jj) = jk      ! Turbocline  
    121120            END DO 
    122121         END DO 
    … …  
    127126            iiki = imld(ji,jj) 
    128127            iikn = nmln(ji,jj) 
    129             imkt = mikt(ji,jj) 
    130             hmld (ji,jj) = ( fsdepw(ji,jj,iiki  ) - fsdepw(ji,jj,imkt ) ) * ssmask(ji,jj)    ! Turbocline depth  
    131             hmlp (ji,jj) = ( fsdepw(ji,jj,iikn  ) - fsdepw(ji,jj,imkt ) ) * ssmask(ji,jj)    ! Mixed layer depth 
    132             hmlpt(ji,jj) = ( fsdept(ji,jj,iikn-1) - fsdepw(ji,jj,imkt ) ) * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer 
     128            hmld (ji,jj) = fsdepw(ji,jj,iiki  ) * ssmask(ji,jj)    ! Turbocline depth  
     129            hmlp (ji,jj) = fsdepw(ji,jj,iikn  ) * ssmask(ji,jj)    ! Mixed layer depth 
     130            hmlpt(ji,jj) = fsdept(ji,jj,iikn-1) * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer 
    133131         END DO 
    134132      END DO 
    135       IF( .NOT.lk_offline ) THEN            ! no need to output in offline mode 
    136          CALL iom_put( "mldr10_1", hmlp )   ! mixed layer depth 
    137          CALL iom_put( "mldkz5"  , hmld )   ! turbocline depth 
     133      ! no need to output in offline mode 
     134      IF( .NOT.lk_offline ) THEN    
     135         IF ( iom_use("mldr10_1") ) THEN 
     136            IF( ln_isfcav ) THEN 
     137               CALL iom_put( "mldr10_1", hmlp - risfdep)   ! mixed layer thickness 
     138            ELSE 
     139               CALL iom_put( "mldr10_1", hmlp )            ! mixed layer depth 
     140            END IF 
     141         END IF 
     142         IF ( iom_use("mldkz5") ) THEN 
     143            IF( ln_isfcav ) THEN 
     144               CALL iom_put( "mldkz5"  , hmld - risfdep )   ! turbocline thickness 
     145            ELSE 
     146               CALL iom_put( "mldkz5"  , hmld )             ! turbocline depth 
     147            END IF 
     148         END IF 
    138149      ENDIF 
    139150       

  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/step.F90

    r6237 r6436  
    338338      ! 
    339339      IF( lrst_oce         )   CALL rst_write( kstp )       ! write output ocean restart file 
     340      IF( ln_sto_eos       )   CALL sto_rst_write( kstp )   ! write restart file for stochastic parameters 
    340341 
    341342#if defined key_agrif 

  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90

    r6333 r6436  
    133133         ! 
    134134         CALL p4z_bio( kt, jnt )   ! Biology 
    135          CALL p4z_sed( kt, jnt )   ! Sedimentation 
    136135         CALL p4z_lys( kt, jnt )   ! Compute CaCO3 saturation 
     136         CALL p4z_sed( kt, jnt )   ! Surface and Bottom boundary conditions 
    137137         CALL p4z_flx( kt, jnt )   ! Compute surface fluxes 
    138138         ! 
Note: See TracChangeset for help on using the changeset viewer.