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 5621 for branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Ignore:
Timestamp:
2015-07-21T13:25:36+02:00 (9 years ago)
Author:
mathiot
Message:

UKMO_ISF: upgrade to NEMO_3.6_STABLE (r5554)

Location:
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/LIM_SRC_3
Files:
22 edited

Legend:

Unmodified
Added
Removed

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r5146 r5621  
    373373   INTEGER          , PUBLIC ::   nlay_s          !: number of snow layers  
    374374   CHARACTER(len=32), PUBLIC ::   cn_icerst_in    !: suffix of ice restart name (input) 
     375   CHARACTER(len=256), PUBLIC ::   cn_icerst_indir !: ice restart input directory 
    375376   CHARACTER(len=32), PUBLIC ::   cn_icerst_out   !: suffix of ice restart name (output) 
     377   CHARACTER(len=256), PUBLIC ::   cn_icerst_outdir!: ice restart output directory 
    376378   LOGICAL          , PUBLIC ::   ln_limdyn       !: flag for ice dynamics (T) or not (F) 
    377379   LOGICAL          , PUBLIC ::   ln_icectl       !: flag for sea-ice points output (T) or not (F) 
    … …  
    394396   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_smv  !: transport of salt content 
    395397   ! 
    396    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_heat_dhc !: snw/ice heat content variation   [W/m2]  
     398   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_heat     !: snw/ice heat content variation   [W/m2]  
     399   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_smvi     !: ice salt content variation   []  
     400   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_vice     !: ice volume variation   [m/s]  
     401   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_vsnw     !: snw volume variation   [m/s]  
    397402   ! 
    398403   !!---------------------------------------------------------------------- 
    … …  
    455460      ALLOCATE( t_s(jpi,jpj,nlay_s,jpl) , e_s(jpi,jpj,nlay_s,jpl) , STAT=ierr(ii) ) 
    456461      ii = ii + 1 
    457       ALLOCATE( t_i(jpi,jpj,nlay_i+1,jpl) , e_i(jpi,jpj,nlay_i+1,jpl) , s_i(jpi,jpj,nlay_i+1,jpl) , STAT=ierr(ii) ) 
     462      ALLOCATE( t_i(jpi,jpj,nlay_i,jpl) , e_i(jpi,jpj,nlay_i,jpl) , s_i(jpi,jpj,nlay_i,jpl) , STAT=ierr(ii) ) 
    458463 
    459464      ! * Moments for advection 
    … …  
    471476         &      STAT=ierr(ii) ) 
    472477      ii = ii + 1 
    473       ALLOCATE( sxe (jpi,jpj,nlay_i+1,jpl) , sye (jpi,jpj,nlay_i+1,jpl) , sxxe(jpi,jpj,nlay_i+1,jpl) ,     & 
    474          &      syye(jpi,jpj,nlay_i+1,jpl) , sxye(jpi,jpj,nlay_i+1,jpl)                           , STAT=ierr(ii) ) 
     478      ALLOCATE( sxe (jpi,jpj,nlay_i,jpl) , sye (jpi,jpj,nlay_i,jpl) , sxxe(jpi,jpj,nlay_i,jpl) ,     & 
     479         &      syye(jpi,jpj,nlay_i,jpl) , sxye(jpi,jpj,nlay_i,jpl)                            , STAT=ierr(ii) ) 
    475480 
    476481      ! * Old values of global variables 
    477482      ii = ii + 1 
    478483      ALLOCATE( v_s_b  (jpi,jpj,jpl) , v_i_b  (jpi,jpj,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) ,     & 
    479          &      a_i_b  (jpi,jpj,jpl) , smv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i+1 ,jpl) ,  & 
    480          &      oa_i_b (jpi,jpj,jpl) , u_ice_b(jpi,jpj)     , v_ice_b(jpi,jpj)             , STAT=ierr(ii) ) 
     484         &      a_i_b  (jpi,jpj,jpl) , smv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) ,     & 
     485         &      oa_i_b (jpi,jpj,jpl) , u_ice_b(jpi,jpj)     , v_ice_b(jpi,jpj)          , STAT=ierr(ii) ) 
    481486       
    482487      ! * Ice thickness distribution variables 
    … …  
    486491      ! * Ice diagnostics 
    487492      ii = ii + 1 
    488       ALLOCATE( diag_trp_vi(jpi,jpj), diag_trp_vs (jpi,jpj), diag_trp_ei  (jpi,jpj),   &  
    489          &      diag_trp_es(jpi,jpj), diag_trp_smv(jpi,jpj), diag_heat_dhc(jpi,jpj),  STAT=ierr(ii) ) 
     493      ALLOCATE( diag_trp_vi(jpi,jpj), diag_trp_vs (jpi,jpj), diag_trp_ei(jpi,jpj),   &  
     494         &      diag_trp_es(jpi,jpj), diag_trp_smv(jpi,jpj), diag_heat  (jpi,jpj),   & 
     495         &      diag_smvi  (jpi,jpj), diag_vice   (jpi,jpj), diag_vsnw  (jpi,jpj), STAT=ierr(ii) ) 
    490496 
    491497      ice_alloc = MAXVAL( ierr(:) ) 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90

    r5123 r5621  
    207207 
    208208      !-- Lateral boundary conditions 
    209       CALL lbc_lnk( psm , 'T',  1. )   ;   CALL lbc_lnk( ps0 , 'T',  1. ) 
    210       CALL lbc_lnk( psx , 'T', -1. )   ;   CALL lbc_lnk( psy , 'T', -1. )      ! caution gradient ==> the sign changes 
    211       CALL lbc_lnk( psxx, 'T',  1. )   ;   CALL lbc_lnk( psyy, 'T',  1. ) 
    212       CALL lbc_lnk( psxy, 'T',  1. ) 
     209      CALL lbc_lnk_multi( psm , 'T',  1., ps0 , 'T',  1.   & 
     210         &              , psx , 'T', -1., psy , 'T', -1.   &   ! caution gradient ==> the sign changes 
     211         &              , psxx, 'T',  1., psyy, 'T',  1.   & 
     212         &              , psxy, 'T',  1. ) 
    213213 
    214214      IF(ln_ctl) THEN 
    … …  
    393393 
    394394      !-- Lateral boundary conditions 
    395       CALL lbc_lnk( psm , 'T',  1. )   ;   CALL lbc_lnk( ps0 , 'T',  1. ) 
    396       CALL lbc_lnk( psx , 'T', -1. )   ;   CALL lbc_lnk( psy , 'T', -1. )      ! caution gradient ==> the sign changes 
    397       CALL lbc_lnk( psxx, 'T',  1. )   ;   CALL lbc_lnk( psyy, 'T',  1. ) 
    398       CALL lbc_lnk( psxy, 'T',  1. ) 
     395      CALL lbc_lnk_multi( psm , 'T',  1.,  ps0 , 'T',  1.   & 
     396         &              , psx , 'T', -1.,  psy , 'T', -1.   &   ! caution gradient ==> the sign changes 
     397         &              , psxx, 'T',  1.,  psyy, 'T',  1.   & 
     398         &              , psxy, 'T',  1. ) 
    399399 
    400400      IF(ln_ctl) THEN 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90

    r5123 r5621  
    88   !!            3.5  ! 2011-02  (G. Madec)  add mpp considerations 
    99   !!             -   ! 2014-05  (C. Rousset) add lim_cons_hsm 
     10   !!             -   ! 2015-03  (C. Rousset) add lim_cons_final 
    1011   !!---------------------------------------------------------------------- 
    1112#if defined key_lim3 
    … …  
    2223   USE lib_mpp        ! MPP library 
    2324   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     25   USE sbc_oce , ONLY : sfx  ! Surface boundary condition: ocean fields 
    2426 
    2527   IMPLICIT NONE 
    … …  
    3032   PUBLIC   lim_cons_check 
    3133   PUBLIC   lim_cons_hsm 
     34   PUBLIC   lim_cons_final 
    3235 
    3336   !!---------------------------------------------------------------------- 
    … …  
    7275      !! ** Method  : Arithmetics 
    7376      !!--------------------------------------------------------------------- 
    74       INTEGER                                  , INTENT(in   ) ::   ksum   !: number of categories 
    75       INTEGER                                  , INTENT(in   ) ::   klay   !: number of vertical layers 
    76       REAL(wp), DIMENSION(jpi,jpj,nlay_i+1,jpl), INTENT(in   ) ::   pin   !: input field 
    77       REAL(wp), DIMENSION(jpi,jpj)             , INTENT(  out) ::   pout   !: output field 
     77      INTEGER                                , INTENT(in   ) ::   ksum   !: number of categories 
     78      INTEGER                                , INTENT(in   ) ::   klay   !: number of vertical layers 
     79      REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl), INTENT(in   ) ::   pin    !: input field 
     80      REAL(wp), DIMENSION(jpi,jpj)           , INTENT(  out) ::   pout   !: output field 
    7881      ! 
    7982      INTEGER ::   jk, jl   ! dummy loop indices 
    … …  
    155158 
    156159   SUBROUTINE lim_cons_hsm( icount, cd_routine, zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b ) 
    157       !!------------------------------------------------------------------- 
    158       !!               ***  ROUTINE lim_cons_hsm *** 
    159       !! 
    160       !! ** Purpose : Test the conservation of heat, salt and mass for each routine 
    161       !! 
    162       !! ** Method  : 
    163       !!--------------------------------------------------------------------- 
    164       INTEGER         , INTENT(in)    :: icount      ! determine wether this is the beggining of the routine (0) or the end (1) 
    165       CHARACTER(len=*), INTENT(in)    :: cd_routine  ! name of the routine 
     160      !!-------------------------------------------------------------------------------------------------------- 
     161      !!                                        ***  ROUTINE lim_cons_hsm *** 
     162      !! 
     163      !! ** Purpose : Test the conservation of heat, salt and mass for each ice routine 
     164      !!                     + test if ice concentration and volume are > 0 
     165      !! 
     166      !! ** Method  : This is an online diagnostics which can be activated with ln_limdiahsb=true 
     167      !!              It prints in ocean.output if there is a violation of conservation at each time-step 
     168      !!              The thresholds (zv_sill, zs_sill, zh_sill) which determine violations are set to 
     169      !!              a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years. 
     170      !!              For salt and heat thresholds, ice is considered to have a salinity of 10  
     171      !!              and a heat content of 3e5 J/kg (=latent heat of fusion)  
     172      !!-------------------------------------------------------------------------------------------------------- 
     173      INTEGER         , INTENT(in)    :: icount        ! determine wether this is the beggining of the routine (0) or the end (1) 
     174      CHARACTER(len=*), INTENT(in)    :: cd_routine    ! name of the routine 
    166175      REAL(wp)        , INTENT(inout) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    167176      REAL(wp)                        :: zvi,   zsmv,   zei,   zfs,   zfw,   zft 
    168177      REAL(wp)                        :: zvmin, zamin, zamax  
    169       REAL(wp)                        :: zconv 
    170  
    171       zconv = 1.e-9 
     178      REAL(wp)                        :: zvtrp, zetrp 
     179      REAL(wp)                        :: zarea, zv_sill, zs_sill, zh_sill 
     180      REAL(wp), PARAMETER             :: zconv = 1.e-9 ! convert W to GW and kg to Mt 
    172181 
    173182      IF( icount == 0 ) THEN 
    174183 
     184         ! salt flux 
    175185         zfs_b  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
    176186            &                  sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:)                                  & 
    177             &                ) *  e12t(:,:) * tmask(:,:,1) ) 
    178  
     187            &                ) *  e12t(:,:) * tmask(:,:,1) * zconv ) 
     188 
     189         ! water flux 
    179190         zfw_b  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  & 
    180191            &                  wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:)    & 
    181             &                ) *  e12t(:,:) * tmask(:,:,1) ) 
    182  
     192            &                ) *  e12t(:,:) * tmask(:,:,1) * zconv ) 
     193 
     194         ! heat flux 
    183195         zft_b  = glob_sum(  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  &  
    184196            &                - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   & 
    185197            &                ) *  e12t(:,:) * tmask(:,:,1) * zconv ) 
    186198 
    187          zvi_b  = glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * e12t(:,:) * tmask(:,:,1) ) 
    188  
    189          zsmv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) ) 
     199         zvi_b  = glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e12t * tmask(:,:,1) * zconv ) 
     200 
     201         zsmv_b = glob_sum( SUM( smv_i * rhoic            , dim=3 ) * e12t * tmask(:,:,1) * zconv ) 
    190202 
    191203         zei_b  = glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) +  & 
    192204            &                 SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )    & 
    193                             ) * e12t(:,:) * tmask(:,:,1) * zconv ) 
     205                            ) * e12t * tmask(:,:,1) * zconv ) 
    194206 
    195207      ELSEIF( icount == 1 ) THEN 
    196208 
     209         ! salt flux 
    197210         zfs  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
    198211            &                sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:)                                  &  
    199             &              ) * e12t(:,:) * tmask(:,:,1) ) - zfs_b 
    200  
     212            &              ) * e12t(:,:) * tmask(:,:,1) * zconv ) - zfs_b 
     213 
     214         ! water flux 
    201215         zfw  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  & 
    202216            &                wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:)    & 
    203             &              ) * e12t(:,:) * tmask(:,:,1) ) - zfw_b 
    204  
     217            &              ) * e12t(:,:) * tmask(:,:,1) * zconv ) - zfw_b 
     218 
     219         ! heat flux 
    205220         zft  = glob_sum(  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  &  
    206221            &              - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   & 
    207222            &              ) * e12t(:,:) * tmask(:,:,1) * zconv ) - zft_b 
    208223  
    209          zvi  = ( glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 )  & 
    210             &                    * e12t(:,:) * tmask(:,:,1) ) - zvi_b ) * r1_rdtice - zfw  
    211  
    212          zsmv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) ) - zsmv_b ) * r1_rdtice + ( zfs * r1_rhoic ) 
     224         ! outputs 
     225         zvi  = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 )  & 
     226            &                    * e12t * tmask(:,:,1) * zconv ) - zvi_b ) * r1_rdtice - zfw ) * rday 
     227 
     228         zsmv = ( ( glob_sum( SUM( smv_i * rhoic            , dim=3 )  & 
     229            &                    * e12t * tmask(:,:,1) * zconv ) - zsmv_b ) * r1_rdtice + zfs ) * rday 
    213230 
    214231         zei  =   glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) +  & 
    215232            &                 SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )    & 
    216             &                ) * e12t(:,:) * tmask(:,:,1) * zconv ) * r1_rdtice - zei_b * r1_rdtice + zft 
     233            &                ) * e12t * tmask(:,:,1) * zconv ) * r1_rdtice - zei_b * r1_rdtice + zft 
     234 
     235         ! zvtrp and zetrp must be close to 0 if the advection scheme is conservative 
     236         zvtrp = glob_sum( ( diag_trp_vi * rhoic + diag_trp_vs * rhosn ) * e12t * tmask(:,:,1) * zconv ) * rday  
     237         zetrp = glob_sum( ( diag_trp_ei         + diag_trp_es         ) * e12t * tmask(:,:,1) * zconv ) 
    217238 
    218239         zvmin = glob_min( v_i ) 
    219240         zamax = glob_max( SUM( a_i, dim=3 ) ) 
    220241         zamin = glob_min( a_i ) 
    221         
     242 
     243         ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice)  
     244         zarea   = glob_sum( SUM( a_i + epsi10, dim=3 ) * e12t * zconv ) ! in 1.e9 m2 
     245         zv_sill = zarea * 2.5e-5 
     246         zs_sill = zarea * 25.e-5 
     247         zh_sill = zarea * 10.e-5 
     248 
    222249         IF(lwp) THEN 
    223             IF ( ABS( zvi    ) >  1.e-4 ) WRITE(numout,*) 'violation volume [kg/day]     (',cd_routine,') = ',(zvi * rday) 
    224             IF ( ABS( zsmv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (',cd_routine,') = ',(zsmv * rday) 
    225             IF ( ABS( zei    ) >  1.e-4 ) WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',(zei) 
    226             IF ( zvmin <  -epsi10       ) WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',(zvmin) 
    227             IF( cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' .AND. zamax > rn_amax+epsi10 ) THEN 
    228                                           WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax 
     250            IF ( ABS( zvi  ) > zv_sill ) WRITE(numout,*) 'violation volume [Mt/day]     (',cd_routine,') = ',zvi 
     251            IF ( ABS( zsmv ) > zs_sill ) WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zsmv 
     252            IF ( ABS( zei  ) > zh_sill ) WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',zei 
     253            IF ( ABS(zvtrp ) > zv_sill .AND. cd_routine == 'limtrp' ) THEN 
     254                                         WRITE(numout,*) 'violation vtrp [Mt/day]       (',cd_routine,') = ',zvtrp 
     255                                         WRITE(numout,*) 'violation etrp [GW]           (',cd_routine,') = ',zetrp 
    229256            ENDIF 
    230             IF ( zamin <  -epsi10       ) WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin 
     257            IF (     zvmin   < -epsi10 ) WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',zvmin 
     258            IF (     zamax   > rn_amax+epsi10 .AND. cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' ) THEN 
     259                                         WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax 
     260            ENDIF 
     261            IF (      zamin  < -epsi10 ) WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin 
    231262         ENDIF 
    232263 
    … …  
    234265 
    235266   END SUBROUTINE lim_cons_hsm 
     267 
     268   SUBROUTINE lim_cons_final( cd_routine ) 
     269      !!--------------------------------------------------------------------------------------------------------- 
     270      !!                                   ***  ROUTINE lim_cons_final *** 
     271      !! 
     272      !! ** Purpose : Test the conservation of heat, salt and mass at the end of each ice time-step 
     273      !! 
     274      !! ** Method  : This is an online diagnostics which can be activated with ln_limdiahsb=true 
     275      !!              It prints in ocean.output if there is a violation of conservation at each time-step 
     276      !!              The thresholds (zv_sill, zs_sill, zh_sill) which determine the violation are set to 
     277      !!              a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years. 
     278      !!              For salt and heat thresholds, ice is considered to have a salinity of 10  
     279      !!              and a heat content of 3e5 J/kg (=latent heat of fusion)  
     280      !!-------------------------------------------------------------------------------------------------------- 
     281      CHARACTER(len=*), INTENT(in)    :: cd_routine    ! name of the routine 
     282      REAL(wp)                        :: zhfx, zsfx, zvfx 
     283      REAL(wp)                        :: zarea, zv_sill, zs_sill, zh_sill 
     284      REAL(wp), PARAMETER             :: zconv = 1.e-9 ! convert W to GW and kg to Mt 
     285 
     286#if ! defined key_bdy 
     287      ! heat flux 
     288      zhfx  = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es - hfx_sub ) * e12t * tmask(:,:,1) * zconv )  
     289      ! salt flux 
     290      zsfx  = glob_sum( ( sfx + diag_smvi ) * e12t * tmask(:,:,1) * zconv ) * rday 
     291      ! water flux 
     292      zvfx  = glob_sum( ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e12t * tmask(:,:,1) * zconv ) * rday 
     293 
     294      ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice)  
     295      zarea   = glob_sum( SUM( a_i + epsi10, dim=3 ) * e12t * zconv ) ! in 1.e9 m2 
     296      zv_sill = zarea * 2.5e-5 
     297      zs_sill = zarea * 25.e-5 
     298      zh_sill = zarea * 10.e-5 
     299 
     300      IF( ABS( zvfx ) > zv_sill ) WRITE(numout,*) 'violation vfx    [Mt/day]       (',cd_routine,')  = ',(zvfx) 
     301      IF( ABS( zsfx ) > zs_sill ) WRITE(numout,*) 'violation sfx    [psu*Mt/day]   (',cd_routine,')  = ',(zsfx) 
     302      IF( ABS( zhfx ) > zh_sill ) WRITE(numout,*) 'violation hfx    [GW]           (',cd_routine,')  = ',(zhfx) 
     303#endif 
     304 
     305   END SUBROUTINE lim_cons_final 
    236306 
    237307#else 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90

    r5125 r5621  
    419419               WRITE(numout,*) ' hfx_in       : ', hfx_in(ji,jj) 
    420420               WRITE(numout,*) ' hfx_out      : ', hfx_out(ji,jj) 
    421                WRITE(numout,*) ' dhc          : ', diag_heat_dhc(ji,jj)               
     421               WRITE(numout,*) ' dhc          : ', diag_heat(ji,jj)               
    422422               WRITE(numout,*) 
    423423               WRITE(numout,*) ' hfx_dyn      : ', hfx_dyn(ji,jj) 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90

    • Property svn:keywords set to Id
    r5123 r5621  
    4040   !!---------------------------------------------------------------------- 
    4141   !! NEMO/OPA 3.4 , NEMO Consortium (2012) 
    42    !! $Id: limdiahsb.F90 3294 2012-10-18 16:44:18Z rblod $ 
     42   !! $Id$ 
    4343   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4444   !!---------------------------------------------------------------------- 
    … …  
    115115      zbg_ihc      = glob_sum( et_i(:,:) * e12t(:,:) * 1.e-20 ) ! ice heat content  [1.e20 J] 
    116116      zbg_shc      = glob_sum( et_s(:,:) * e12t(:,:) * 1.e-20 ) ! snow heat content [1.e20 J] 
    117       zbg_hfx_dhc  = glob_sum( diag_heat_dhc(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
     117      zbg_hfx_dhc  = glob_sum( diag_heat(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    118118      zbg_hfx_spr  = glob_sum( hfx_spr(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    119119 
    … …  
    245245         WRITE(numout,*) '~~~~~~~~~~~~' 
    246246      ENDIF 
    247  
    248       ! ---------------------------------- ! 
    249       ! 2 - initial conservation variables ! 
    250       ! ---------------------------------- ! 
    251       !frc_vol = 0._wp                                          ! volume       trend due to forcing 
    252       !frc_sal = 0._wp                                          ! salt content   -    -   -    -          
    253       !bg_grme = 0._wp                                          ! ice growth + melt volume trend 
    254247      ! 
    255248      CALL lim_diahsb_rst( nstart, 'READ' )  !* read or initialize all required files 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90

    r5123 r5621  
    1313   !!---------------------------------------------------------------------- 
    1414   !!   lim_hdf       : diffusion trend on sea-ice variable 
     15   !!   lim_hdf_init  : initialisation of diffusion trend on sea-ice variable 
    1516   !!---------------------------------------------------------------------- 
    1617   USE dom_oce        ! ocean domain 
    … …  
    2627   PRIVATE 
    2728 
    28    PUBLIC   lim_hdf     ! called by lim_trp 
     29   PUBLIC   lim_hdf         ! called by lim_trp 
     30   PUBLIC   lim_hdf_init    ! called by sbc_lim_init 
    2931 
    3032   LOGICAL  ::   linit = .TRUE.                             ! initialization flag (set to flase after the 1st call) 
     33   INTEGER  ::   nn_convfrq                                 !:  convergence check frequency of the Crant-Nicholson scheme 
    3134   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   efact   ! metric coefficient 
    3235 
    … …  
    121124         CALL lbc_lnk( zrlx, 'T', 1. )                   ! lateral boundary condition 
    122125         ! 
    123          zconv = 0._wp                                   ! convergence test 
    124          DO jj = 2, jpjm1 
    125             DO ji = fs_2, fs_jpim1 
    126                zconv = MAX( zconv, ABS( zrlx(ji,jj) - ptab(ji,jj) )  ) 
    127             END DO 
    128          END DO 
    129          IF( lk_mpp )   CALL mpp_max( zconv )            ! max over the global domain 
     126         IF ( MOD( iter, nn_convfrq ) == 0 )  THEN    ! convergence test every nn_convfrq iterations (perf. optimization) 
     127            zconv = 0._wp 
     128            DO jj = 2, jpjm1 
     129               DO ji = fs_2, fs_jpim1 
     130                  zconv = MAX( zconv, ABS( zrlx(ji,jj) - ptab(ji,jj) )  ) 
     131               END DO 
     132            END DO 
     133            IF( lk_mpp )   CALL mpp_max( zconv )      ! max over the global domain 
     134         ENDIF 
    130135         ! 
    131136         ptab(:,:) = zrlx(:,:) 
    … …  
    162167   END SUBROUTINE lim_hdf 
    163168 
     169    
     170   SUBROUTINE lim_hdf_init 
     171      !!------------------------------------------------------------------- 
     172      !!                  ***  ROUTINE lim_hdf_init  *** 
     173      !! 
     174      !! ** Purpose : Initialisation of horizontal diffusion of sea-ice  
     175      !! 
     176      !! ** Method  : Read the namicehdf namelist 
     177      !! 
     178      !! ** input   : Namelist namicehdf 
     179      !!------------------------------------------------------------------- 
     180      INTEGER  ::   ios                 ! Local integer output status for namelist read 
     181      NAMELIST/namicehdf/ nn_convfrq 
     182      !!------------------------------------------------------------------- 
     183      ! 
     184      IF(lwp) THEN 
     185         WRITE(numout,*) 
     186         WRITE(numout,*) 'lim_hdf : Ice horizontal diffusion' 
     187         WRITE(numout,*) '~~~~~~~' 
     188      ENDIF 
     189      ! 
     190      REWIND( numnam_ice_ref )              ! Namelist namicehdf in reference namelist : Ice horizontal diffusion 
     191      READ  ( numnam_ice_ref, namicehdf, IOSTAT = ios, ERR = 901) 
     192901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicehdf in reference namelist', lwp ) 
     193 
     194      REWIND( numnam_ice_cfg )              ! Namelist namicehdf in configuration namelist : Ice horizontal diffusion 
     195      READ  ( numnam_ice_cfg, namicehdf, IOSTAT = ios, ERR = 902 ) 
     196902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicehdf in configuration namelist', lwp ) 
     197      IF(lwm) WRITE ( numoni, namicehdf ) 
     198      ! 
     199      IF(lwp) THEN                          ! control print 
     200         WRITE(numout,*) 
     201         WRITE(numout,*)'   Namelist of ice parameters for ice horizontal diffusion computation ' 
     202         WRITE(numout,*)'      convergence check frequency of the Crant-Nicholson scheme   nn_convfrq   = ', nn_convfrq 
     203      ENDIF 
     204      ! 
     205   END SUBROUTINE lim_hdf_init 
    164206#else 
    165207   !!---------------------------------------------------------------------- 
    166208   !!   Default option          Dummy module           NO LIM sea-ice model 
    167209   !!---------------------------------------------------------------------- 
    168 CONTAINS 
    169    SUBROUTINE lim_hdf         ! Empty routine 
    170    END SUBROUTINE lim_hdf 
    171210#endif 
    172211 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90

    r5123 r5621  
    117117 
    118118      ! basal temperature (considered at freezing point) 
    119       t_bo(:,:) = ( eos_fzp( tsn(:,:,1,jp_sal) ) + rt0 ) * tmask(:,:,1)  
     119      CALL eos_fzp( sss_m(:,:), t_bo(:,:) ) 
     120      t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1)  
    120121 
    121122      IF( ln_iceini ) THEN 
    … …  
    127128      DO jj = 1, jpj                                       ! ice if sst <= t-freez + ttest 
    128129         DO ji = 1, jpi 
    129             IF( ( tsn(ji,jj,1,jp_tem)  - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst ) THEN  
     130            IF( ( sst_m(ji,jj)  - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst ) THEN  
    130131               zswitch(ji,jj) = 0._wp * tmask(ji,jj,1)    ! no ice 
    131132            ELSE                                                                                    
    … …  
    314315            DO ji = 1, jpi 
    315316               a_i(ji,jj,jl)   = zswitch(ji,jj) * za_i_ini (jl,zhemis(ji,jj))  ! concentration 
    316                ht_i(ji,jj,jl)  = zswitch(ji,jj) * zh_i_ini(jl,zhemis(ji,jj))  ! ice thickness 
     317               ht_i(ji,jj,jl)  = zswitch(ji,jj) * zh_i_ini(jl,zhemis(ji,jj))   ! ice thickness 
    317318               ht_s(ji,jj,jl)  = ht_i(ji,jj,jl) * ( zht_s_ini( zhemis(ji,jj) ) / zht_i_ini( zhemis(ji,jj) ) )  ! snow depth 
    318                sm_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) !+ ( 1._wp - zswitch(ji,jj) ) * rn_simin ! salinity 
    319                o_i(ji,jj,jl)   = zswitch(ji,jj) * 1._wp + ( 1._wp - zswitch(ji,jj) ) ! age 
     319               sm_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj))     ! salinity 
     320               o_i(ji,jj,jl)   = zswitch(ji,jj) * 1._wp                        ! age (1 day) 
    320321               t_su(ji,jj,jl)  = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp 
    321322 
    … …  
    333334               smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content 
    334335               oa_i(ji,jj,jl)  = o_i(ji,jj,jl) * a_i(ji,jj,jl)               ! age content 
    335             END DO ! ji 
    336          END DO ! jj 
    337       END DO ! jl 
     336            END DO 
     337         END DO 
     338      END DO 
    338339 
    339340      ! Snow temperature and heat content 
    … …  
    348349                   ! Mutliply by volume, and divide by number of layers to get heat content in J/m2 
    349350                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s 
    350                END DO ! ji 
    351             END DO ! jj 
    352          END DO ! jl 
    353       END DO ! jk 
     351               END DO 
     352            END DO 
     353         END DO 
     354      END DO 
    354355 
    355356      ! Ice salinity, temperature and heat content 
    … …  
    369370                   ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 
    370371                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i 
    371                END DO ! ji 
    372             END DO ! jj 
    373          END DO ! jl 
    374       END DO ! jk 
     372               END DO 
     373            END DO 
     374         END DO 
     375      END DO 
    375376 
    376377      tn_ice (:,:,:) = t_su (:,:,:) 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90

    r5134 r5621  
    127127      REAL(wp) ::   za, zfac              ! local scalar 
    128128      CHARACTER (len = 15) ::   fieldid 
    129       REAL(wp), POINTER, DIMENSION(:,:) ::   closing_net     ! net rate at which area is removed    (1/s) 
    130                                                              ! (ridging ice area - area of new ridges) / dt 
    131       REAL(wp), POINTER, DIMENSION(:,:) ::   divu_adv        ! divu as implied by transport scheme  (1/s) 
    132       REAL(wp), POINTER, DIMENSION(:,:) ::   opning          ! rate of opening due to divergence/shear 
    133       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 
     129      REAL(wp), POINTER, DIMENSION(:,:)   ::   closing_net     ! net rate at which area is removed    (1/s) 
     130                                                               ! (ridging ice area - area of new ridges) / dt 
     131      REAL(wp), POINTER, DIMENSION(:,:)   ::   divu_adv        ! divu as implied by transport scheme  (1/s) 
     132      REAL(wp), POINTER, DIMENSION(:,:)   ::   opning          ! rate of opening due to divergence/shear 
     133      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 
    137137      ! 
    138138      INTEGER, PARAMETER ::   nitermax = 20     
    … …  
    142142      IF( nn_timing == 1 )  CALL timing_start('limitd_me') 
    143143 
    144       CALL wrk_alloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 
     144      CALL wrk_alloc( jpi,jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 
    145145 
    146146      IF(ln_ctl) THEN 
    … …  
    153153      ! conservation test 
    154154      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     155 
     156      CALL lim_var_zapsmall 
     157      CALL lim_var_glo2eqv            ! equivalent variables, requested for rafting 
    155158 
    156159      !-----------------------------------------------------------------------------! 
    … …  
    235238               ! Reduce the closing rate if more than 100% of the open water  
    236239               ! would be removed.  Reduce the opening rate proportionately. 
    237                IF ( ato_i(ji,jj) > epsi10 .AND. athorn(ji,jj,0) > 0.0 ) THEN 
    238                   za = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice 
    239                   IF ( za > ato_i(ji,jj)) THEN 
    240                      zfac = ato_i(ji,jj) / za 
    241                      closing_gross(ji,jj) = closing_gross(ji,jj) * zfac 
    242                      opning(ji,jj) = opning(ji,jj) * zfac 
    243                   ENDIF 
     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 
    244245               ENDIF 
    245246 
    … …  
    251252         ! Reduce the closing rate if more than 100% of any ice category  
    252253         ! would be removed.  Reduce the opening rate proportionately. 
    253  
    254254         DO jl = 1, jpl 
    255255            DO jj = 1, jpj 
    256256               DO ji = 1, jpi 
    257                   IF ( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp )THEN 
    258                      za = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice 
    259                      IF ( za  >  a_i(ji,jj,jl) ) THEN 
    260                         zfac = a_i(ji,jj,jl) / za 
    261                         closing_gross(ji,jj) = closing_gross(ji,jj) * zfac 
    262                         opning       (ji,jj) = opning       (ji,jj) * zfac 
    263                      ENDIF 
     257                  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 ) 
     260                     closing_gross(ji,jj) = closing_gross(ji,jj) * zfac 
     261                     opning       (ji,jj) = opning       (ji,jj) * zfac 
    264262                  ENDIF 
    265263               END DO 
    … …  
    368366      ENDIF 
    369367 
    370       ! updates 
    371       CALL lim_var_glo2eqv 
    372       CALL lim_var_zapsmall 
    373368      CALL lim_var_agg( 1 )  
    374369 
    … …  
    377372      !-----------------------------------------------------------------------------! 
    378373      IF(ln_ctl) THEN  
     374         CALL lim_var_glo2eqv 
     375 
    379376         CALL prt_ctl_info(' ') 
    380377         CALL prt_ctl_info(' - Cell values : ') 
    … …  
    531528         DO jj = 2, jpjm1 
    532529            DO ji = 2, jpim1 
    533                IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > epsi10) THEN ! ice is present 
     530               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp) THEN  
    534531                  zworka(ji,jj) = ( 4.0 * strength(ji,jj)              & 
    535532                     &                  + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) &   
    … …  
    566563         DO jj = 1, jpj - 1 
    567564            DO ji = 1, jpi - 1 
    568                IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > epsi10) THEN       ! ice is present 
     565               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp) THEN  
    569566                  numts_rm = 1 ! number of time steps for the running mean 
    570567                  IF ( strp1(ji,jj) > 0.0 ) numts_rm = numts_rm + 1 
    … …  
    637634 
    638635      Gsum(:,:,-1) = 0._wp 
    639  
    640       DO jj = 1, jpj 
    641          DO ji = 1, jpi 
    642             IF( ato_i(ji,jj) > epsi10 ) THEN   ;   Gsum(ji,jj,0) = ato_i(ji,jj) 
    643             ELSE                               ;   Gsum(ji,jj,0) = 0._wp 
    644             ENDIF 
    645          END DO 
    646       END DO 
     636      Gsum(:,:,0 ) = ato_i(:,:) 
    647637 
    648638      ! for each value of h, you have to add ice concentration then 
    649639      DO jl = 1, jpl 
    650          DO jj = 1, jpj  
    651             DO ji = 1, jpi 
    652                IF( a_i(ji,jj,jl) > epsi10 ) THEN   ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl) 
    653                ELSE                                ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) 
    654                ENDIF 
    655             END DO 
    656          END DO 
     640         Gsum(:,:,jl) = Gsum(:,:,jl-1) + a_i(:,:,jl) 
    657641      END DO 
    658642 
    … …  
    828812      LOGICAL, PARAMETER ::   l_conservation_check = .true.  ! if true, check conservation (useful for debugging) 
    829813      ! 
    830       LOGICAL ::   neg_ato_i      ! flag for ato_i(i,j) < -puny 
    831       LOGICAL ::   large_afrac    ! flag for afrac > 1 
    832       LOGICAL ::   large_afrft    ! flag for afrac > 1 
    833814      INTEGER ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices 
    834815      INTEGER ::   ij                ! horizontal index, combines i and j loops 
    … …  
    850831      REAL(wp), POINTER, DIMENSION(:,:) ::   ardg1 , ardg2    ! area of ice ridged & new ridges 
    851832      REAL(wp), POINTER, DIMENSION(:,:) ::   vsrdg , esrdg    ! snow volume & energy of ridging ice 
    852       REAL(wp), POINTER, DIMENSION(:,:) ::   oirdg1, oirdg2   ! areal age content of ridged & rifging ice 
    853833      REAL(wp), POINTER, DIMENSION(:,:) ::   dhr   , dhr2     ! hrmax - hrmin  &  hrmax^2 - hrmin^2 
    854834 
    … …  
    859839      REAL(wp), POINTER, DIMENSION(:,:) ::   srdg2   ! sal*volume of new ridges 
    860840      REAL(wp), POINTER, DIMENSION(:,:) ::   smsw    ! sal*volume of water trapped into ridges 
     841      REAL(wp), POINTER, DIMENSION(:,:) ::   oirdg1, oirdg2   ! ice age of ice ridged 
    861842 
    862843      REAL(wp), POINTER, DIMENSION(:,:) ::   afrft            ! fraction of category area rafted 
    … …  
    864845      REAL(wp), POINTER, DIMENSION(:,:) ::   virft , vsrft    ! ice & snow volume of rafting ice 
    865846      REAL(wp), POINTER, DIMENSION(:,:) ::   esrft , smrft    ! snow energy & salinity of rafting ice 
    866       REAL(wp), POINTER, DIMENSION(:,:) ::   oirft1, oirft2   ! areal age content of rafted ice & rafting ice 
     847      REAL(wp), POINTER, DIMENSION(:,:) ::   oirft1, oirft2   ! ice age of ice rafted 
    867848 
    868849      REAL(wp), POINTER, DIMENSION(:,:,:) ::   eirft      ! ice energy of rafting ice 
    … …  
    872853      !!---------------------------------------------------------------------- 
    873854 
    874       CALL wrk_alloc( (jpi+1)*(jpj+1),      indxi, indxj ) 
    875       CALL wrk_alloc( jpi, jpj,             vice_init, vice_final, eice_init, eice_final ) 
    876       CALL wrk_alloc( jpi, jpj,             afrac, fvol , ardg1, ardg2, vsrdg, esrdg, oirdg1, oirdg2, dhr, dhr2 ) 
    877       CALL wrk_alloc( jpi, jpj,             vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw ) 
    878       CALL wrk_alloc( jpi, jpj,             afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 ) 
    879       CALL wrk_alloc( jpi, jpj, jpl,        aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init ) 
    880       CALL wrk_alloc( jpi, jpj, nlay_i+1,      eirft, erdg1, erdg2, ersw ) 
    881       CALL wrk_alloc( jpi, jpj, nlay_i+1, jpl, eicen_init ) 
     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 ) 
    882863 
    883864      ! Conservation check 
    … …  
    898879      ! 1) Compute change in open water area due to closing and opening. 
    899880      !------------------------------------------------------------------------------- 
    900  
    901       neg_ato_i = .false. 
    902  
    903881      DO jj = 1, jpj 
    904882         DO ji = 1, jpi 
    905883            ato_i(ji,jj) = ato_i(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice        & 
    906884               &                        + opning(ji,jj)                          * rdt_ice 
    907             IF( ato_i(ji,jj) < -epsi10 ) THEN 
    908                neg_ato_i = .TRUE. 
    909             ELSEIF( ato_i(ji,jj) < 0._wp ) THEN    ! roundoff error 
     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 
    910888               ato_i(ji,jj) = 0._wp 
    911889            ENDIF 
    912890         END DO 
    913891      END DO 
    914  
    915       ! if negative open water area alert it 
    916       IF( neg_ato_i .AND. lwp ) THEN       ! there is a bug 
    917          DO jj = 1, jpj  
    918             DO ji = 1, jpi 
    919                IF( ato_i(ji,jj) < -epsi10 ) THEN  
    920                   WRITE(numout,*) ''   
    921                   WRITE(numout,*) 'Ridging error: ato_i < 0' 
    922                   WRITE(numout,*) 'ato_i : ', ato_i(ji,jj) 
    923                ENDIF 
    924             END DO 
    925          END DO 
    926       ENDIF 
    927892 
    928893      !----------------------------------------------------------------- 
    929894      ! 2) Save initial state variables 
    930895      !----------------------------------------------------------------- 
    931  
    932       DO jl = 1, jpl 
    933          aicen_init(:,:,jl) = a_i(:,:,jl) 
    934          vicen_init(:,:,jl) = v_i(:,:,jl) 
    935          vsnwn_init(:,:,jl) = v_s(:,:,jl) 
    936          ! 
    937          smv_i_init(:,:,jl) = smv_i(:,:,jl) 
    938          oa_i_init (:,:,jl) = oa_i (:,:,jl) 
    939       END DO 
    940  
    941       esnwn_init(:,:,:) = e_s(:,:,1,:) 
    942  
    943       DO jl = 1, jpl   
    944          DO jk = 1, nlay_i 
    945             eicen_init(:,:,jk,jl) = e_i(:,:,jk,jl) 
    946          END DO 
    947       END DO 
     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 (:,:,:) 
    948903 
    949904      ! 
    … …  
    972927         END DO 
    973928 
    974          large_afrac = .false. 
    975          large_afrft = .false. 
    976  
    977929         DO ij = 1, icells 
    978930            ji = indxi(ij) 
    … …  
    988940            arft2(ji,jj) = arft1(ji,jj) / kraft 
    989941 
    990             oirdg1(ji,jj)= aridge(ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice 
    991             oirft1(ji,jj)= araft (ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice 
    992             oirdg2(ji,jj)= oirdg1(ji,jj) / krdg(ji,jj,jl1) 
    993             oirft2(ji,jj)= oirft1(ji,jj) / kraft 
    994  
    995942            !--------------------------------------------------------------- 
    996943            ! 3.3) Compute ridging /rafting fractions, make sure afrac <=1  
    … …  
    1000947            afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting 
    1001948 
    1002             IF (afrac(ji,jj) > kamax + epsi10) THEN  !riging 
    1003                large_afrac = .true. 
    1004             ELSEIF (afrac(ji,jj) > kamax) THEN  ! roundoff error 
     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 
    1005952               afrac(ji,jj) = kamax 
    1006953            ENDIF 
    1007             IF (afrft(ji,jj) > kamax + epsi10) THEN !rafting 
    1008                large_afrft = .true. 
    1009             ELSEIF (afrft(ji,jj) > kamax) THEN  ! roundoff error 
     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 
    1010958               afrft(ji,jj) = kamax 
    1011959            ENDIF 
    … …  
    1019967            vsw  (ji,jj) = vrdg1(ji,jj) * rn_por_rdg 
    1020968 
    1021             vsrdg(ji,jj) = vsnwn_init(ji,jj,jl1) * afrac(ji,jj) 
    1022             esrdg(ji,jj) = esnwn_init(ji,jj,jl1) * afrac(ji,jj) 
    1023             srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) 
    1024             srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) !! MV HC 2014 this line seems useless 
     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)  
    1025974 
    1026975            ! rafting volumes, heat contents ... 
    1027             virft(ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj) 
    1028             vsrft(ji,jj) = vsnwn_init(ji,jj,jl1) * afrft(ji,jj) 
    1029             esrft(ji,jj) = esnwn_init(ji,jj,jl1) * afrft(ji,jj) 
    1030             smrft(ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj)  
     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  
    1031982 
    1032983            ! substract everything 
    1033             a_i(ji,jj,jl1)   = a_i(ji,jj,jl1)   - ardg1(ji,jj)  - arft1(ji,jj) 
    1034             v_i(ji,jj,jl1)   = v_i(ji,jj,jl1)   - vrdg1(ji,jj)  - virft(ji,jj) 
    1035             v_s(ji,jj,jl1)   = v_s(ji,jj,jl1)   - vsrdg(ji,jj)  - vsrft(ji,jj) 
    1036             e_s(ji,jj,1,jl1) = e_s(ji,jj,1,jl1) - esrdg(ji,jj)  - esrft(ji,jj) 
     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) 
    1037989            oa_i(ji,jj,jl1)  = oa_i(ji,jj,jl1)  - oirdg1(ji,jj) - oirft1(ji,jj) 
    1038             smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - srdg1(ji,jj)  - smrft(ji,jj) 
    1039990 
    1040991            !----------------------------------------------------------------- 
    1041992            ! 3.5) Compute properties of new ridges 
    1042993            !----------------------------------------------------------------- 
    1043             !------------- 
     994            !--------- 
    1044995            ! Salinity 
    1045             !------------- 
     996            !--------- 
    1046997            smsw(ji,jj)  = vsw(ji,jj) * sss_m(ji,jj)                      ! salt content of seawater frozen in voids !! MV HC2014 
    1047998            srdg2(ji,jj) = srdg1(ji,jj) + smsw(ji,jj)                     ! salt content of new ridge 
    … …  
    10501001             
    10511002            sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ji,jj) * rhoic * r1_rdtice 
    1052             wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ji,jj) * rhoic * r1_rdtice   ! gurvan: increase in ice volume du to seawater frozen in voids              
     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              
    10531004 
    10541005            !------------------------------------             
    … …  
    11341085         ENDIF 
    11351086 
    1136          IF( large_afrac .AND. lwp ) THEN   ! there is a bug 
    1137             DO ij = 1, icells 
    1138                ji = indxi(ij) 
    1139                jj = indxj(ij) 
    1140                IF( afrac(ji,jj) > kamax + epsi10 ) THEN  
    1141                   WRITE(numout,*) '' 
    1142                   WRITE(numout,*) ' ardg > a_i' 
    1143                   WRITE(numout,*) ' ardg, aicen_init : ', ardg1(ji,jj), aicen_init(ji,jj,jl1) 
    1144                ENDIF 
    1145             END DO 
    1146          ENDIF 
    1147          IF( large_afrft .AND. lwp ) THEN  ! there is a bug 
    1148             DO ij = 1, icells 
    1149                ji = indxi(ij) 
    1150                jj = indxj(ij) 
    1151                IF( afrft(ji,jj) > kamax + epsi10 ) THEN  
    1152                   WRITE(numout,*) '' 
    1153                   WRITE(numout,*) ' arft > a_i' 
    1154                   WRITE(numout,*) ' arft, aicen_init : ', arft1(ji,jj), aicen_init(ji,jj,jl1) 
    1155                ENDIF 
    1156             END DO 
    1157          ENDIF 
    1158  
    11591087         !------------------------------------------------------------------------------- 
    11601088         ! 4) Add area, volume, and energy of new ridge to each category jl2 
    … …  
    11901118               oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirdg2(ji,jj) * farea 
    11911119 
    1192             END DO ! ij 
     1120            END DO 
    11931121 
    11941122            ! Transfer ice energy to category jl2 by ridging 
    … …  
    12171145                  e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + esrft (ji,jj) * rn_fsnowrft 
    12181146                  smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + smrft (ji,jj)     
    1219                   oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirft2(ji,jj)     
     1147                  oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirft2(ji,jj) 
    12201148               ENDIF 
    12211149               ! 
    … …  
    12571185      ENDIF 
    12581186      ! 
    1259       CALL wrk_dealloc( (jpi+1)*(jpj+1),      indxi, indxj ) 
    1260       CALL wrk_dealloc( jpi, jpj,             vice_init, vice_final, eice_init, eice_final ) 
    1261       CALL wrk_dealloc( jpi, jpj,             afrac, fvol , ardg1, ardg2, vsrdg, esrdg, oirdg1, oirdg2, dhr, dhr2 ) 
    1262       CALL wrk_dealloc( jpi, jpj,             vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw ) 
    1263       CALL wrk_dealloc( jpi, jpj,             afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 ) 
    1264       CALL wrk_dealloc( jpi, jpj, jpl,        aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init ) 
    1265       CALL wrk_dealloc( jpi, jpj, nlay_i+1,      eirft, erdg1, erdg2, ersw ) 
    1266       CALL wrk_dealloc( jpi, jpj, nlay_i+1, jpl, eicen_init ) 
     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 ) 
    12671195      ! 
    12681196   END SUBROUTINE lim_itd_me_ridgeshift 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90

    r5134 r5621  
    9191      !!------------------------------------------------------------------ 
    9292 
    93       CALL wrk_alloc( jpi,jpj, zremap_flag )    ! integer 
    94       CALL wrk_alloc( jpi,jpj,jpl-1, zdonor )   ! integer 
     93      CALL wrk_alloc( jpi,jpj, zremap_flag ) 
     94      CALL wrk_alloc( jpi,jpj,jpl-1, zdonor ) 
    9595      CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 
    9696      CALL wrk_alloc( jpi,jpj,jpl-1, zdaice, zdvice )    
    9797      CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )    
    9898      CALL wrk_alloc( (jpi+1)*(jpj+1), zvetamin, zvetamax )    
    99       CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j )   ! integer  
     99      CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j )  
    100100      CALL wrk_alloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final ) 
    101101 
    … …  
    128128               rswitch           = MAX( 0.0, SIGN( 1.0, a_i(ji,jj,jl) - epsi10 ) )     !0 if no ice and 1 if yes 
    129129               ht_i(ji,jj,jl)    = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * rswitch 
    130                rswitch           = MAX( 0.0, SIGN( 1.0, a_i_b(ji,jj,jl) - epsi10) )    !0 if no ice and 1 if yes 
     130               rswitch           = MAX( 0.0, SIGN( 1.0, a_i_b(ji,jj,jl) - epsi10) ) 
    131131               zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * rswitch 
    132 !clem               IF( a_i(ji,jj,jl) > epsi10 )   zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl)  
    133                zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl)  
     132               IF( a_i(ji,jj,jl) > epsi10 )   zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl) ! clem: useless IF statement?  
    134133            END DO 
    135134         END DO 
    … …  
    173172            ! 
    174173            zhbnew(ii,ij,jl) = hi_max(jl) 
    175             IF ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN 
     174            IF    ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN 
    176175               !interpolate between adjacent category growth rates 
    177176               zslope           = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) ) 
    178177               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_b(ii,ij,jl) ) 
    179             ELSEIF ( a_i_b(ii,ij,jl) > epsi10) THEN 
     178            ELSEIF( a_i_b(ii,ij,jl) > epsi10) THEN 
    180179               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) 
    181             ELSEIF ( a_i_b(ii,ij,jl+1) > epsi10) THEN 
     180            ELSEIF( a_i_b(ii,ij,jl+1) > epsi10) THEN 
    182181               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 
    183182            ENDIF 
    … …  
    188187            ii = nind_i(ji) 
    189188            ij = nind_j(ji) 
    190             IF( a_i(ii,ij,jl) > epsi10 .AND. ht_i(ii,ij,jl) >= zhbnew(ii,ij,jl) ) THEN 
     189 
     190            ! clem: we do not want ht_i to be too close to either HR or HL otherwise a division by nearly 0 is possible  
     191            ! in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 
     192            IF    ( a_i(ii,ij,jl  ) > epsi10 .AND. ht_i(ii,ij,jl  ) > ( zhbnew(ii,ij,jl) - epsi10 ) ) THEN 
    191193               zremap_flag(ii,ij) = 0 
    192             ELSEIF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) <= zhbnew(ii,ij,jl) ) THEN 
     194            ELSEIF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) < ( zhbnew(ii,ij,jl) + epsi10 ) ) THEN 
    193195               zremap_flag(ii,ij) = 0 
    194196            ENDIF 
    195197 
    196198            !- 4.3 Check that each zhbnew does not exceed maximal values hi_max   
     199            IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 
    197200            IF( zhbnew(ii,ij,jl) > hi_max(jl+1) ) zremap_flag(ii,ij) = 0 
    198             IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 
     201            ! clem bug: why is not the following instead? 
     202            !!IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 
     203            !!IF( zhbnew(ii,ij,jl) > hi_max(jl  ) ) zremap_flag(ii,ij) = 0 
     204  
    199205         END DO 
    200206 
    … …  
    220226      DO jj = 1, jpj 
    221227         DO ji = 1, jpi 
    222             zhb0(ji,jj) = hi_max(0) ! 0eme 
    223             zhb1(ji,jj) = hi_max(1) ! 1er 
    224  
    225             zhbnew(ji,jj,klbnd-1) = 0._wp 
     228            zhb0(ji,jj) = hi_max(0) 
     229            zhb1(ji,jj) = hi_max(1) 
    226230 
    227231            IF( a_i(ji,jj,kubnd) > epsi10 ) THEN 
    228232               zhbnew(ji,jj,kubnd) = MAX( hi_max(kubnd-1), 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1) ) 
    229233            ELSE 
    230                zhbnew(ji,jj,kubnd) = hi_max(kubnd)   
    231                !!? clem bug: since hi_max(jpl)=99, this limit is very high  
    232                !!? but I think it is erased in fitline subroutine  
     234!clem bug               zhbnew(ji,jj,kubnd) = hi_max(kubnd)   
     235               zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) ! not used anyway 
     236            ENDIF 
     237 
     238            ! clem: we do not want ht_i_b to be too close to either HR or HL otherwise a division by nearly 0 is possible  
     239            ! in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 
     240            IF    ( zht_i_b(ji,jj,klbnd) < ( zhb0(ji,jj) + epsi10 ) )  THEN 
     241               zremap_flag(ji,jj) = 0 
     242            ELSEIF( zht_i_b(ji,jj,klbnd) > ( zhb1(ji,jj) - epsi10 ) )  THEN 
     243               zremap_flag(ji,jj) = 0 
    233244            ENDIF 
    234245 
    … …  
    249260 
    250261         IF( a_i(ii,ij,klbnd) > epsi10 ) THEN 
     262 
    251263            zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category 
    252             IF( zdh0 < 0.0 ) THEN !remove area from category 1 
     264 
     265            IF( zdh0 < 0.0 ) THEN      !remove area from category 1 
    253266               zdh0 = MIN( -zdh0, hi_max(klbnd) ) 
    254  
    255267               !Integrate g(1) from 0 to dh0 to estimate area melted 
    256268               zetamax = MIN( zdh0, hR(ii,ij,klbnd) ) - hL(ii,ij,klbnd) 
     269 
    257270               IF( zetamax > 0.0 ) THEN 
    258                   zx1  = zetamax 
    259                   zx2  = 0.5 * zetamax * zetamax  
    260                   zda0 = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1 !ice area removed 
    261                   ! Constrain new thickness <= ht_i 
    262                   zdamax = a_i(ii,ij,klbnd) * (1.0 - ht_i(ii,ij,klbnd) / zht_i_b(ii,ij,klbnd) ) ! zdamax > 0 
    263                   !ice area lost due to melting of thin ice 
    264                   zda0   = MIN( zda0, zdamax ) 
    265  
     271                  zx1    = zetamax 
     272                  zx2    = 0.5 * zetamax * zetamax  
     273                  zda0   = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1                        ! ice area removed 
     274                  zdamax = a_i(ii,ij,klbnd) * (1.0 - ht_i(ii,ij,klbnd) / zht_i_b(ii,ij,klbnd) ) ! Constrain new thickness <= ht_i                 
     275                  zda0   = MIN( zda0, zdamax )                                                  ! ice area lost due to melting  
     276                                                                                                !     of thin ice (zdamax > 0) 
    266277                  ! Remove area, conserving volume 
    267278                  ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 ) 
    … …  
    270281               ENDIF 
    271282 
    272             ELSE ! if ice accretion ! a_i > epsi10; zdh0 > 0 
     283            ELSE ! if ice accretion zdh0 > 0 
     284               ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1) 
    273285               zhbnew(ii,ij,klbnd-1) = MIN( zdh0, hi_max(klbnd) )  
    274                ! zhbnew was 0, and is shifted to the right to account for thin ice 
    275                ! growth in openwater (F0 = f1) 
    276             ENDIF ! zdh0  
    277  
    278          ENDIF ! a_i > epsi10 
     286            ENDIF 
     287 
     288         ENDIF 
    279289 
    280290      END DO 
    … …  
    304314 
    305315            IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! transfer from jl to jl+1 
    306  
    307316               ! left and right integration limits in eta space 
    308317               zvetamin(ji) = MAX( hi_max(jl), hL(ii,ij,jl) ) - hL(ii,ij,jl) 
    309                zvetamax(ji) = MIN (zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl) 
     318               zvetamax(ji) = MIN( zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl) 
    310319               zdonor(ii,ij,jl) = jl 
    311320 
    312             ELSE  ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl 
    313  
     321            ELSE                                    ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl 
    314322               ! left and right integration limits in eta space 
    315323               zvetamin(ji) = 0.0 
    … …  
    317325               zdonor(ii,ij,jl) = jl + 1 
    318326 
    319             ENDIF  ! zhbnew(jl) > hi_max(jl) 
     327            ENDIF 
    320328 
    321329            zetamax = MAX( zvetamax(ji), zvetamin(ji) ) ! no transfer if etamax < etamin 
    … …  
    334342 
    335343         END DO 
    336       END DO ! jl klbnd -> kubnd - 1 
     344      END DO 
    337345 
    338346      !!---------------------------------------------------------------------------------------------- 
    … …  
    376384      ENDIF 
    377385 
    378       CALL wrk_dealloc( jpi,jpj, zremap_flag )    ! integer 
    379       CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor )   ! integer 
     386      CALL wrk_dealloc( jpi,jpj, zremap_flag ) 
     387      CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor ) 
    380388      CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 
    381389      CALL wrk_dealloc( jpi,jpj,jpl-1, zdaice, zdvice )    
    382390      CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )    
    383391      CALL wrk_dealloc( (jpi+1)*(jpj+1), zvetamin, zvetamax )    
    384       CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j )   ! integer  
     392      CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j )  
    385393      CALL wrk_dealloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final ) 
    386394 
    … …  
    407415      INTEGER , DIMENSION(jpi,jpj), INTENT(in   ) ::   zremap_flag  ! 
    408416      ! 
    409       INTEGER ::   ji,jj           ! horizontal indices 
     417      INTEGER  ::   ji,jj        ! horizontal indices 
    410418      REAL(wp) ::   zh13         ! HbL + 1/3 * (HbR - HbL) 
    411419      REAL(wp) ::   zh23         ! HbL + 2/3 * (HbR - HbL) 
    … …  
    414422      !!------------------------------------------------------------------ 
    415423      ! 
    416       ! 
    417424      DO jj = 1, jpj 
    418425         DO ji = 1, jpi 
    419426            ! 
    420427            IF( zremap_flag(ji,jj) == 1 .AND. a_i(ji,jj,num_cat) > epsi10   & 
    421                &                        .AND. hice(ji,jj)        > 0._wp     ) THEN 
     428               &                        .AND. hice(ji,jj)        > 0._wp ) THEN 
    422429 
    423430               ! Initialize hL and hR 
    424  
    425431               hL(ji,jj) = HbL(ji,jj) 
    426432               hR(ji,jj) = HbR(ji,jj) 
    427433 
    428434               ! Change hL or hR if hice falls outside central third of range 
    429  
    430435               zh13 = 1.0 / 3.0 * ( 2.0 * hL(ji,jj) + hR(ji,jj) ) 
    431436               zh23 = 1.0 / 3.0 * ( hL(ji,jj) + 2.0 * hR(ji,jj) ) 
    … …  
    436441 
    437442               ! Compute coefficients of g(eta) = g0 + g1*eta 
    438  
    439443               zdhr = 1._wp / (hR(ji,jj) - hL(ji,jj)) 
    440444               zwk1 = 6._wp * a_i(ji,jj,num_cat) * zdhr 
    … …  
    443447               g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 ) 
    444448               ! 
    445             ELSE                   ! remap_flag = .false. or a_i < epsi10  
     449            ELSE  ! remap_flag = .false. or a_i < epsi10  
    446450               hL(ji,jj) = 0._wp 
    447451               hR(ji,jj) = 0._wp 
    448452               g0(ji,jj) = 0._wp 
    449453               g1(ji,jj) = 0._wp 
    450             ENDIF                  ! a_i > epsi10 
     454            ENDIF 
    451455            ! 
    452456         END DO 
    … …  
    472476 
    473477      INTEGER ::   ji, jj, jl, jl2, jl1, jk   ! dummy loop indices 
    474       INTEGER ::   ii, ij          ! indices when changing from 2D-1D is done 
     478      INTEGER ::   ii, ij                     ! indices when changing from 2D-1D is done 
    475479 
    476480      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zaTsfn 
    … …  
    485489      INTEGER, POINTER, DIMENSION(:) ::   nind_i, nind_j   ! compressed indices for i/j directions 
    486490 
    487       INTEGER ::   nbrem             ! number of cells with ice to transfer 
    488  
    489       LOGICAL ::   zdaice_negative         ! true if daice < -puny 
    490       LOGICAL ::   zdvice_negative         ! true if dvice < -puny 
    491       LOGICAL ::   zdaice_greater_aicen    ! true if daice > aicen 
    492       LOGICAL ::   zdvice_greater_vicen    ! true if dvice > vicen 
     491      INTEGER  ::   nbrem             ! number of cells with ice to transfer 
    493492      !!------------------------------------------------------------------ 
    494493 
    495494      CALL wrk_alloc( jpi,jpj,jpl, zaTsfn ) 
    496495      CALL wrk_alloc( jpi,jpj, zworka ) 
    497       CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j )   ! integer 
     496      CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 
    498497 
    499498      !---------------------------------------------------------------------------------------------- 
    … …  
    505504      END DO 
    506505 
    507 !clem: I think the following is wrong (if enabled, it creates negative concentration/volume around -epsi10) 
    508 !      !---------------------------------------------------------------------------------------------- 
    509 !      ! 2) Check for daice or dvice out of range, allowing for roundoff error 
    510 !      !---------------------------------------------------------------------------------------------- 
    511 !      ! Note: zdaice < 0 or zdvice < 0 usually happens when category jl 
    512 !      ! has a small area, with h(n) very close to a boundary.  Then 
    513 !      ! the coefficients of g(h) are large, and the computed daice and 
    514 !      ! dvice can be in error. If this happens, it is best to transfer 
    515 !      ! either the entire category or nothing at all, depending on which 
    516 !      ! side of the boundary hice(n) lies. 
    517 !      !----------------------------------------------------------------- 
    518 !      DO jl = klbnd, kubnd-1 
    519 ! 
    520 !         zdaice_negative = .false. 
    521 !         zdvice_negative = .false. 
    522 !         zdaice_greater_aicen = .false. 
    523 !         zdvice_greater_vicen = .false. 
    524 ! 
    525 !         DO jj = 1, jpj 
    526 !            DO ji = 1, jpi 
    527 ! 
    528 !               IF (zdonor(ji,jj,jl) > 0) THEN 
    529 !                  jl1 = zdonor(ji,jj,jl) 
    530 ! 
    531 !                  IF (zdaice(ji,jj,jl) < 0.0) THEN 
    532 !                     IF (zdaice(ji,jj,jl) > -epsi10) THEN 
    533 !                        IF ( ( jl1 == jl   .AND. ht_i(ji,jj,jl1) >  hi_max(jl) ) .OR.  & 
    534 !                             ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN 
    535 !                           zdaice(ji,jj,jl) = a_i(ji,jj,jl1)  ! shift entire category 
    536 !                           zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 
    537 !                        ELSE 
    538 !                           zdaice(ji,jj,jl) = 0.0 ! shift no ice 
    539 !                           zdvice(ji,jj,jl) = 0.0 
    540 !                        ENDIF 
    541 !                     ELSE 
    542 !                        zdaice_negative = .true. 
    543 !                     ENDIF 
    544 !                  ENDIF 
    545 ! 
    546 !                  IF (zdvice(ji,jj,jl) < 0.0) THEN 
    547 !                     IF (zdvice(ji,jj,jl) > -epsi10 ) THEN 
    548 !                        IF ( ( jl1 == jl   .AND. ht_i(ji,jj,jl1) >  hi_max(jl) ) .OR.  & 
    549 !                             ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN 
    550 !                           zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 
    551 !                           zdvice(ji,jj,jl) = v_i(ji,jj,jl1)  
    552 !                        ELSE 
    553 !                           zdaice(ji,jj,jl) = 0.0    ! shift no ice 
    554 !                           zdvice(ji,jj,jl) = 0.0 
    555 !                        ENDIF 
    556 !                    ELSE 
    557 !                       zdvice_negative = .true. 
    558 !                    ENDIF 
    559 !                 ENDIF 
    560 ! 
    561 !                  ! If daice is close to aicen, set daice = aicen. 
    562 !                  IF (zdaice(ji,jj,jl) > a_i(ji,jj,jl1) - epsi10 ) THEN 
    563 !                     IF (zdaice(ji,jj,jl) < a_i(ji,jj,jl1)+epsi10) THEN 
    564 !                        zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 
    565 !                        zdvice(ji,jj,jl) = v_i(ji,jj,jl1)  
    566 !                    ELSE 
    567 !                       zdaice_greater_aicen = .true. 
    568 !                    ENDIF 
    569 !                  ENDIF 
    570 ! 
    571 !                  IF (zdvice(ji,jj,jl) > v_i(ji,jj,jl1)-epsi10) THEN 
    572 !                     IF (zdvice(ji,jj,jl) < v_i(ji,jj,jl1)+epsi10) THEN 
    573 !                        zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 
    574 !                        zdvice(ji,jj,jl) = v_i(ji,jj,jl1)  
    575 !                     ELSE 
    576 !                        zdvice_greater_vicen = .true. 
    577 !                     ENDIF 
    578 !                  ENDIF 
    579 ! 
    580 !               ENDIF               ! donor > 0 
    581 !            END DO 
    582 !         END DO 
    583 ! 
    584 !      END DO 
    585 !clem 
    586506      !------------------------------------------------------------------------------- 
    587       ! 3) Transfer volume and energy between categories 
     507      ! 2) Transfer volume and energy between categories 
    588508      !------------------------------------------------------------------------------- 
    589509 
    … …  
    605525 
    606526            jl1 = zdonor(ii,ij,jl) 
    607             rswitch       = MAX( 0._wp , SIGN( 1._wp , v_i(ii,ij,jl1) - epsi20 ) ) 
    608             zworka(ii,ij) = zdvice(ii,ij,jl) / MAX( v_i(ii,ij,jl1), epsi20 ) * rswitch 
     527            rswitch       = MAX( 0._wp , SIGN( 1._wp , v_i(ii,ij,jl1) - epsi10 ) ) 
     528            zworka(ii,ij) = zdvice(ii,ij,jl) / MAX( v_i(ii,ij,jl1), epsi10 ) * rswitch 
    609529            IF( jl1 == jl) THEN   ;   jl2 = jl1+1 
    610530            ELSE                  ;   jl2 = jl  
    … …  
    614534            ! Ice areas 
    615535            !-------------- 
    616  
    617536            a_i(ii,ij,jl1) = a_i(ii,ij,jl1) - zdaice(ii,ij,jl) 
    618537            a_i(ii,ij,jl2) = a_i(ii,ij,jl2) + zdaice(ii,ij,jl) 
    … …  
    621540            ! Ice volumes 
    622541            !-------------- 
    623  
    624542            v_i(ii,ij,jl1) = v_i(ii,ij,jl1) - zdvice(ii,ij,jl)  
    625543            v_i(ii,ij,jl2) = v_i(ii,ij,jl2) + zdvice(ii,ij,jl) 
    … …  
    628546            ! Snow volumes 
    629547            !-------------- 
    630  
    631548            zdvsnow        = v_s(ii,ij,jl1) * zworka(ii,ij) 
    632549            v_s(ii,ij,jl1) = v_s(ii,ij,jl1) - zdvsnow 
    … …  
    636553            ! Snow heat content   
    637554            !-------------------- 
    638  
    639555            zdesnow            = e_s(ii,ij,1,jl1) * zworka(ii,ij) 
    640556            e_s(ii,ij,1,jl1)   = e_s(ii,ij,1,jl1) - zdesnow 
    … …  
    644560            ! Ice age  
    645561            !-------------- 
    646  
    647562            zdo_aice           = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl) 
    648563            oa_i(ii,ij,jl1)    = oa_i(ii,ij,jl1) - zdo_aice 
    … …  
    652567            ! Ice salinity 
    653568            !-------------- 
    654  
    655569            zdsm_vice          = smv_i(ii,ij,jl1) * zworka(ii,ij) 
    656570            smv_i(ii,ij,jl1)   = smv_i(ii,ij,jl1) - zdsm_vice 
    … …  
    660574            ! Surface temperature 
    661575            !--------------------- 
    662  
    663576            zdaTsf             = t_su(ii,ij,jl1) * zdaice(ii,ij,jl) 
    664577            zaTsfn(ii,ij,jl1)  = zaTsfn(ii,ij,jl1) - zdaTsf 
    … …  
    711624      CALL wrk_dealloc( jpi,jpj,jpl, zaTsfn ) 
    712625      CALL wrk_dealloc( jpi,jpj, zworka ) 
    713       CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j )   ! integer 
     626      CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 
    714627      ! 
    715628   END SUBROUTINE lim_itd_shiftice 
    … …  
    737650      REAL(wp), POINTER, DIMENSION(:,:) ::   vt_s_init, vt_s_final   ! snow volume summed over categories 
    738651      !!------------------------------------------------------------------ 
    739       !! clem 2014/04: be carefull, rebining does not conserve salt(maybe?) => the difference is taken into account in limupdate 
    740652       
    741653      CALL wrk_alloc( jpi,jpj,jpl, zdonor )   ! interger 
    … …  
    845757         ENDIF 
    846758 
    847 !         ! clem-change begin: why not doing that? 
    848 !         DO jj = 1, jpj 
    849 !            DO ji = 1, jpi 
    850 !               IF( a_i(ji,jj,jl+1) > epsi10 .AND. ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 
    851 !                  ht_i(ji,jj,jl+1) = hi_max(jl) + epsi10 
    852 !                  a_i (ji,jj,jl+1) = v_i(ji,jj,jl+1) / ht_i(ji,jj,jl+1)  
    853 !               ENDIF 
    854 !            END DO 
    855 !         END DO 
    856          ! clem-change end 
    857  
    858759      END DO 
    859760 
    … …  
    872773      ENDIF 
    873774      ! 
    874       CALL wrk_dealloc( jpi,jpj,jpl, zdonor )   ! interger 
     775      CALL wrk_dealloc( jpi,jpj,jpl, zdonor ) 
    875776      CALL wrk_dealloc( jpi,jpj,jpl, zdaice, zdvice ) 
    876777      CALL wrk_dealloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final ) 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

    r5123 r5621  
    377377            END DO 
    378378         END DO 
    379          CALL lbc_lnk( v_ice1  , 'U', -1. )   ;   CALL lbc_lnk( u_ice2  , 'V', -1. )      ! lateral boundary cond. 
    380  
     379 
     380         CALL lbc_lnk_multi( v_ice1, 'U', -1., u_ice2, 'V', -1. )      ! lateral boundary cond. 
     381          
    381382         DO jj = k_j1+1, k_jpj-1 
    382383            DO ji = fs_2, fs_jpim1 
    … …  
    412413            END DO 
    413414         END DO 
    414          CALL lbc_lnk( zs1 , 'T', 1. )   ;   CALL lbc_lnk( zs2, 'T', 1. ) 
    415          CALL lbc_lnk( zs12, 'F', 1. ) 
    416  
     415 
     416         CALL lbc_lnk_multi( zs1 , 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 
     417  
    417418         ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) 
    418419         DO jj = k_j1+1, k_jpj-1 
    … …  
    570571      END DO 
    571572 
    572       CALL lbc_lnk( u_ice(:,:), 'U', -1. )  
    573       CALL lbc_lnk( v_ice(:,:), 'V', -1. )  
     573      CALL lbc_lnk_multi( u_ice(:,:), 'U', -1., v_ice(:,:), 'V', -1. ) 
     574 
    574575#if defined key_agrif && defined key_lim2 
    575576      CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'U' ) 
    … …  
    595596      END DO 
    596597 
    597       CALL lbc_lnk( u_ice2(:,:), 'V', -1. )  
    598       CALL lbc_lnk( v_ice1(:,:), 'U', -1. ) 
     598      CALL lbc_lnk_multi( u_ice2(:,:), 'V', -1., v_ice1(:,:), 'U', -1. ) 
    599599 
    600600      ! Recompute delta, shear and div, inputs for mechanical redistribution  
    … …  
    643643 
    644644      ! Lateral boundary condition 
    645       CALL lbc_lnk( divu_i (:,:), 'T', 1. ) 
    646       CALL lbc_lnk( delta_i(:,:), 'T', 1. ) 
    647       ! CALL lbc_lnk( shear_i(:,:), 'F', 1. ) 
    648       CALL lbc_lnk( shear_i(:,:), 'T', 1. ) 
     645      CALL lbc_lnk_multi( divu_i (:,:), 'T', 1., delta_i(:,:), 'T', 1.,  shear_i(:,:), 'T', 1. ) 
    649646 
    650647      ! * Store the stress tensor for the next time step 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90

    r5128 r5621  
    5555      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step define as a character 
    5656      CHARACTER(LEN=50)   ::   clname   ! ice output restart file name 
     57      CHARACTER(len=256)  ::   clpath   ! full path to ice output restart file  
    5758      !!---------------------------------------------------------------------- 
    5859      ! 
    … …  
    6465      IF( kt == nitrst - 2*nn_fsbc + 1 .OR. nstock == nn_fsbc    & 
    6566         &                             .OR. ( kt == nitend - nn_fsbc + 1 .AND. .NOT. lrst_ice ) ) THEN 
    66          ! beware of the format used to write kt (default is i8.8, that should be large enough...) 
    67          IF( nitrst > 99999999 ) THEN   ;   WRITE(clkt, *       ) nitrst 
    68          ELSE                           ;   WRITE(clkt, '(i8.8)') nitrst 
     67         IF( nitrst <= nitend .AND. nitrst > 0 ) THEN 
     68            ! beware of the format used to write kt (default is i8.8, that should be large enough...) 
     69            IF( nitrst > 99999999 ) THEN   ;   WRITE(clkt, *       ) nitrst 
     70            ELSE                           ;   WRITE(clkt, '(i8.8)') nitrst 
     71            ENDIF 
     72            ! create the file 
     73            clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_"//TRIM(cn_icerst_out) 
     74            clpath = TRIM(cn_icerst_outdir)  
     75            IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath)//'/' 
     76            IF(lwp) THEN 
     77               WRITE(numout,*) 
     78               SELECT CASE ( jprstlib ) 
     79               CASE ( jprstdimg ) 
     80                  WRITE(numout,*) '             open ice restart binary file: ',TRIM(clpath)//clname 
     81               CASE DEFAULT 
     82                  WRITE(numout,*) '             open ice restart NetCDF file: ',TRIM(clpath)//clname 
     83               END SELECT 
     84               IF( kt == nitrst - 2*nn_fsbc + 1 ) THEN    
     85                  WRITE(numout,*)         '             kt = nitrst - 2*nn_fsbc + 1 = ', kt,' date= ', ndastp 
     86               ELSE   ;   WRITE(numout,*) '             kt = '                         , kt,' date= ', ndastp 
     87               ENDIF 
     88            ENDIF 
     89            ! 
     90            CALL iom_open( TRIM(clpath)//TRIM(clname), numriw, ldwrt = .TRUE., kiolib = jprstlib ) 
     91            lrst_ice = .TRUE. 
    6992         ENDIF 
    70          ! create the file 
    71          clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_"//TRIM(cn_icerst_out) 
    72          IF(lwp) THEN 
    73             WRITE(numout,*) 
    74             SELECT CASE ( jprstlib ) 
    75             CASE ( jprstdimg )   ;   WRITE(numout,*) '             open ice restart binary file: '//clname 
    76             CASE DEFAULT         ;   WRITE(numout,*) '             open ice restart NetCDF file: '//clname 
    77             END SELECT 
    78             IF( kt == nitrst - 2*nn_fsbc + 1 ) THEN    
    79                WRITE(numout,*)         '             kt = nitrst - 2*nn_fsbc + 1 = ', kt,' date= ', ndastp 
    80             ELSE   ;   WRITE(numout,*) '             kt = '                         , kt,' date= ', ndastp 
    81             ENDIF 
    82          ENDIF 
    83          ! 
    84          CALL iom_open( clname, numriw, ldwrt = .TRUE., kiolib = jprstlib ) 
    85          lrst_ice = .TRUE. 
    8693      ENDIF 
    8794      ! 
    … …  
    143150         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) 
    144151      END DO 
    145        
     152 
    146153      DO jl = 1, jpl  
    147154         WRITE(zchar,'(I1)') jl 
    … …  
    327334        ! eventually read netcdf file (monobloc)  for restarting on different number of processors 
    328335        ! if {cn_icerst_in}.nc exists, then set jlibalt to jpnf90 
    329         INQUIRE( FILE = TRIM(cn_icerst_in)//'.nc', EXIST = llok ) 
     336        INQUIRE( FILE = TRIM(cn_icerst_indir)//'/'//TRIM(cn_icerst_in)//'.nc', EXIST = llok ) 
    330337        IF ( llok ) THEN ; jlibalt = jpnf90  ; ELSE ; jlibalt = jprstlib ; ENDIF 
    331338      ENDIF 
    332339 
    333       CALL iom_open ( cn_icerst_in, numrir, kiolib = jprstlib ) 
     340      CALL iom_open ( TRIM(cn_icerst_indir)//'/'//cn_icerst_in, numrir, kiolib = jprstlib ) 
    334341 
    335342      CALL iom_get( numrir, 'nn_fsbc', zfice ) 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90

    r5146 r5621  
    3030   USE sbc_oce          ! Surface boundary condition: ocean fields 
    3131   USE sbccpl 
    32    USE oce       , ONLY : fraqsr_1lev, sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass 
     32   USE oce       , ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass 
    3333   USE albedo           ! albedo parameters 
    3434   USE lbclnk           ! ocean lateral boundary condition - MPP exchanges 
    … …  
    4242   USE domvvl           ! Variable volume 
    4343   USE limctl 
     44   USE limcons 
    4445 
    4546   IMPLICIT NONE 
    … …  
    9394      !!              - fr_i    : ice fraction 
    9495      !!              - tn_ice  : sea-ice surface temperature 
    95       !!              - alb_ice : sea-ice albedo (lk_cpl=T) 
     96      !!              - alb_ice : sea-ice albedo (only useful in coupled mode) 
    9697      !! 
    9798      !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90. 
    … …  
    100101      !!              The ref should be Rousset et al., 2015 
    101102      !!--------------------------------------------------------------------- 
    102       INTEGER, INTENT(in) ::   kt                                   ! number of iteration 
    103       INTEGER  ::   ji, jj, jl, jk                                  ! dummy loop indices 
    104       REAL(wp) ::   zemp                                            ! local scalars 
    105       REAL(wp) ::   zf_mass                                         ! Heat flux associated with mass exchange ice->ocean (W.m-2) 
    106       REAL(wp) ::   zfcm1                                           ! New solar flux received by the ocean 
     103      INTEGER, INTENT(in) ::   kt                                  ! number of iteration 
     104      INTEGER  ::   ji, jj, jl, jk                                 ! dummy loop indices 
     105      REAL(wp) ::   zqmass                                         ! Heat flux associated with mass exchange ice->ocean (W.m-2) 
     106      REAL(wp) ::   zqsr                                           ! New solar flux received by the ocean 
    107107      ! 
    108108      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb_cs, zalb_os     ! 2D/3D workspace 
    … …  
    110110 
    111111      ! make calls for heat fluxes before it is modified 
    112       IF( iom_use('qsr_oce') )   CALL iom_put( "qsr_oce" , qsr(:,:) * pfrld(:,:) )   !     solar flux at ocean surface 
    113       IF( iom_use('qns_oce') )   CALL iom_put( "qns_oce" , qns(:,:) * pfrld(:,:) )   ! non-solar flux at ocean surface 
    114       IF( iom_use('qsr_ice') )   CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) ) !     solar flux at ice surface 
    115       IF( iom_use('qns_ice') )   CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) ) ! non-solar flux at ice surface 
    116       IF( iom_use('qtr_ice') )   CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) ) !     solar flux transmitted thru ice 
    117       IF( iom_use('qt_oce' ) )   CALL iom_put( "qt_oce"  , ( qsr(:,:) + qns(:,:) ) * pfrld(:,:) )   
    118       IF( iom_use('qt_ice' ) )   CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) ) * a_i_b(:,:,:), dim=3 ) ) 
     112      IF( iom_use('qsr_oce') )   CALL iom_put( "qsr_oce" , qsr_oce(:,:) * pfrld(:,:) )                                   !     solar flux at ocean surface 
     113      IF( iom_use('qns_oce') )   CALL iom_put( "qns_oce" , qns_oce(:,:) * pfrld(:,:) + qemp_oce(:,:) )                   ! non-solar flux at ocean surface 
     114      IF( iom_use('qsr_ice') )   CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux at ice surface 
     115      IF( iom_use('qns_ice') )   CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) ! non-solar flux at ice surface 
     116      IF( iom_use('qtr_ice') )   CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux transmitted thru ice 
     117      IF( iom_use('qt_oce' ) )   CALL iom_put( "qt_oce"  , ( qsr_oce(:,:) + qns_oce(:,:) ) * pfrld(:,:) + qemp_oce(:,:) )   
     118      IF( iom_use('qt_ice' ) )   CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) )   & 
     119         &                                                      * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) 
     120      IF( iom_use('qemp_oce' ) )   CALL iom_put( "qemp_oce"  , qemp_oce(:,:) )   
     121      IF( iom_use('qemp_ice' ) )   CALL iom_put( "qemp_ice"  , qemp_ice(:,:) )   
    119122 
    120123      ! pfrld is the lead fraction at the previous time step (actually between TRP and THD) 
    … …  
    125128            !      heat flux at the ocean surface      ! 
    126129            !------------------------------------------! 
    127             ! Solar heat flux reaching the ocean = zfcm1 (W.m-2)  
     130            ! Solar heat flux reaching the ocean = zqsr (W.m-2)  
    128131            !--------------------------------------------------- 
    129             IF( lk_cpl ) THEN  
    130                !!! LIM2 version zqsr = qsr_tot(ji,jj) + ( fstric(ji,jj) - qsr_ice(ji,jj,1) ) * ( 1.0 - pfrld(ji,jj) ) 
    131                zfcm1 = qsr_tot(ji,jj) 
    132                DO jl = 1, jpl 
    133                   zfcm1 = zfcm1 + ( ftr_ice(ji,jj,jl) - qsr_ice(ji,jj,jl) ) * a_i_b(ji,jj,jl) 
    134                END DO 
    135             ELSE 
    136                !!! LIM2 version zqsr = pfrld(ji,jj) * qsr(ji,jj)  + ( 1.  - pfrld(ji,jj) ) * fstric(ji,jj) 
    137                zfcm1   = pfrld(ji,jj) * qsr(ji,jj) 
    138                DO jl = 1, jpl 
    139                   zfcm1   = zfcm1 + a_i_b(ji,jj,jl) * ftr_ice(ji,jj,jl) 
    140                END DO 
    141             ENDIF 
     132            zqsr = qsr_tot(ji,jj) 
     133            DO jl = 1, jpl 
     134               zqsr = zqsr - a_i_b(ji,jj,jl) * (  qsr_ice(ji,jj,jl) - ftr_ice(ji,jj,jl) )  
     135            END DO 
    142136 
    143137            ! Total heat flux reaching the ocean = hfx_out (W.m-2)  
    144138            !--------------------------------------------------- 
    145             zf_mass        = hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) ! heat flux from snow is 0 (T=0 degC) 
    146             hfx_out(ji,jj) = hfx_out(ji,jj) + zf_mass + zfcm1 
     139            zqmass         = hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) ! heat flux from snow is 0 (T=0 degC) 
     140            hfx_out(ji,jj) = hfx_out(ji,jj) + zqmass + zqsr 
    147141 
    148142            ! Add the residual from heat diffusion equation (W.m-2) 
    … …  
    152146            ! New qsr and qns used to compute the oceanic heat flux at the next time step 
    153147            !--------------------------------------------------- 
    154             qsr(ji,jj) = zfcm1                                       
    155             qns(ji,jj) = hfx_out(ji,jj) - zfcm1               
     148            qsr(ji,jj) = zqsr                                       
     149            qns(ji,jj) = hfx_out(ji,jj) - zqsr               
    156150 
    157151            !------------------------------------------! 
    … …  
    166160            !                     Even if i see Ice melting as a FW and SALT flux 
    167161            !         
    168             !  computing freshwater exchanges at the ice/ocean interface 
    169             IF( lk_cpl ) THEN  
    170                 zemp =   emp_tot(ji,jj)                                    &   ! net mass flux over grid cell 
    171                    &   - emp_ice(ji,jj) * ( 1._wp - pfrld(ji,jj) )         &   ! minus the mass flux intercepted by sea ice 
    172                    &   + sprecip(ji,jj) * ( pfrld(ji,jj) - pfrld(ji,jj)**rn_betas )   ! 
    173             ELSE 
    174                zemp =   emp(ji,jj)     *           pfrld(ji,jj)            &   ! evaporation over oceanic fraction 
    175                   &   - tprecip(ji,jj) * ( 1._wp - pfrld(ji,jj) )          &   ! all precipitation reach the ocean 
    176                   &   + sprecip(ji,jj) * ( 1._wp - pfrld(ji,jj)**rn_betas )       ! except solid precip intercepted by sea-ice 
    177             ENDIF 
    178  
    179162            ! mass flux from ice/ocean 
    180163            wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj)   & 
    … …  
    182165 
    183166            ! mass flux at the ocean/ice interface 
    184             fmmflx(ji,jj) = - wfx_ice(ji,jj) * r1_rdtice                    ! F/M mass flux save at least for biogeochemical model 
    185             emp(ji,jj)    = zemp - wfx_ice(ji,jj) - wfx_snw(ji,jj)       ! mass flux + F/M mass flux (always ice/ocean mass exchange) 
     167            fmmflx(ji,jj) = - ( wfx_ice(ji,jj) + wfx_snw(ji,jj) ) * r1_rdtice  ! F/M mass flux save at least for biogeochemical model 
     168            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) 
    186169             
    187170         END DO 
    … …  
    212195      tn_ice(:,:,:) = t_su(:,:,:)           ! Ice surface temperature                       
    213196 
    214       !------------------------------------------------! 
    215       !    Snow/ice albedo (only if sent to coupler)   ! 
    216       !------------------------------------------------! 
    217       IF( lk_cpl ) THEN          ! coupled case 
    218  
    219             CALL wrk_alloc( jpi, jpj, jpl, zalb_cs, zalb_os ) 
    220  
    221             CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os )  ! cloud-sky and overcast-sky ice albedos 
    222  
    223             alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
    224  
    225             CALL wrk_dealloc( jpi, jpj, jpl, zalb_cs, zalb_os ) 
    226  
    227       ENDIF 
    228  
    229       IF( ln_icectl )   CALL lim_prt( kt, iiceprt, jiceprt, 3, ' - Final state lim_sbc - ' )   ! control print 
     197      !------------------------------------------------------------------------! 
     198      !    Snow/ice albedo (only if sent to coupler, useless in forced mode)   ! 
     199      !------------------------------------------------------------------------! 
     200      CALL wrk_alloc( jpi, jpj, jpl, zalb_cs, zalb_os )     
     201      CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os )  ! cloud-sky and overcast-sky ice albedos 
     202      alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
     203      CALL wrk_dealloc( jpi, jpj, jpl, zalb_cs, zalb_os ) 
     204 
     205      ! conservation test 
     206      IF( ln_limdiahsb ) CALL lim_cons_final( 'limsbc' ) 
     207 
     208      ! control prints 
     209      IF( ln_icectl )   CALL lim_prt( kt, iiceprt, jiceprt, 3, ' - Final state lim_sbc - ' ) 
    230210 
    231211      IF(ln_ctl) THEN 
    … …  
    341321            sice_0(:,:) = 2._wp 
    342322         END WHERE 
    343       ENDIF 
    344        
    345       IF( .NOT. ln_rstart ) THEN 
    346          fraqsr_1lev(:,:) = 1._wp 
    347323      ENDIF 
    348324      ! 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r5146 r5621  
    2222   USE phycst         ! physical constants 
    2323   USE dom_oce        ! ocean space and time domain variables 
    24    USE oce     , ONLY : fraqsr_1lev  
    2524   USE ice            ! LIM: sea-ice variables 
    2625   USE sbc_oce        ! Surface boundary condition: ocean fields 
    … …  
    2827   USE thd_ice        ! LIM thermodynamic sea-ice variables 
    2928   USE dom_ice        ! LIM sea-ice domain 
    30    USE domvvl         ! domain: variable volume level 
    3129   USE limthd_dif     ! LIM: thermodynamics, vertical diffusion 
    3230   USE limthd_dh      ! LIM: thermodynamics, ice and snow thickness variation 
    … …  
    5048   PRIVATE 
    5149 
    52    PUBLIC   lim_thd        ! called by limstp module 
    53    PUBLIC   lim_thd_init   ! called by sbc_lim_init 
     50   PUBLIC   lim_thd         ! called by limstp module 
     51   PUBLIC   lim_thd_init    ! called by sbc_lim_init 
    5452 
    5553   !! * Substitutions 
    … …  
    8987      REAL(wp) :: zfric_u, zqld, zqfr 
    9088      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    91       REAL(wp), PARAMETER :: zfric_umin = 0._wp        ! lower bound for the friction velocity (cice value=5.e-04) 
    92       REAL(wp), PARAMETER :: zch        = 0.0057_wp    ! heat transfer coefficient 
    93       ! 
    94       REAL(wp), POINTER, DIMENSION(:,:) ::  zqsr, zqns 
     89      REAL(wp), PARAMETER :: zfric_umin = 0._wp           ! lower bound for the friction velocity (cice value=5.e-04) 
     90      REAL(wp), PARAMETER :: zch        = 0.0057_wp       ! heat transfer coefficient 
     91      ! 
    9592      !!------------------------------------------------------------------- 
    96       CALL wrk_alloc( jpi, jpj, zqsr, zqns ) 
    9793 
    9894      IF( nn_timing == 1 )  CALL timing_start('limthd') 
    … …  
    10197      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    10298 
     99      CALL lim_var_glo2eqv 
    103100      !------------------------------------------------------------------------! 
    104101      ! 1) Initialization of some variables                                    ! 
    … …  
    135132      ! 2) Partial computation of forcing for the thermodynamic sea ice model.      ! 
    136133      !-----------------------------------------------------------------------------! 
    137  
    138       !--- Ocean solar and non solar fluxes to be used in zqld 
    139       IF ( .NOT. lk_cpl ) THEN   ! --- forced case, fluxes to the lead are the same as over the ocean 
    140          ! 
    141          zqsr(:,:) = qsr(:,:)      ; zqns(:,:) = qns(:,:) 
    142          ! 
    143       ELSE                       ! --- coupled case, fluxes to the lead are total - intercepted 
    144          ! 
    145          zqsr(:,:) = qsr_tot(:,:)  ; zqns(:,:) = qns_tot(:,:) 
    146          ! 
    147          DO jl = 1, jpl 
    148             DO jj = 1, jpj 
    149                DO ji = 1, jpi 
    150                   zqsr(ji,jj) = zqsr(ji,jj) - qsr_ice(ji,jj,jl) * a_i_b(ji,jj,jl) 
    151                   zqns(ji,jj) = zqns(ji,jj) - qns_ice(ji,jj,jl) * a_i_b(ji,jj,jl) 
    152                END DO 
    153             END DO 
    154          END DO 
    155          ! 
    156       ENDIF 
    157  
    158134      DO jj = 1, jpj 
    159135         DO ji = 1, jpi 
    … …  
    166142            !           !  temperature and turbulent mixing (McPhee, 1992) 
    167143            ! 
    168  
    169144            ! --- Energy received in the lead, zqld is defined everywhere (J.m-2) --- ! 
    170             ! REMARK valid at least in forced mode from clem 
    171             ! precip is included in qns but not in qns_ice 
    172             IF ( lk_cpl ) THEN 
    173                zqld =  tmask(ji,jj,1) * rdt_ice *  & 
    174                   &    (   zqsr(ji,jj) * fraqsr_1lev(ji,jj) + zqns(ji,jj)               &   ! pfrld already included in coupled mode 
    175                   &    + ( pfrld(ji,jj)**rn_betas - pfrld(ji,jj) ) * sprecip(ji,jj)  *     &   ! heat content of precip 
    176                   &      ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus )   & 
    177                   &    + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 ) ) 
    178             ELSE 
    179                zqld =  tmask(ji,jj,1) * rdt_ice *  & 
    180                   &      ( pfrld(ji,jj) * ( zqsr(ji,jj) * fraqsr_1lev(ji,jj) + zqns(ji,jj) )    & 
    181                   &    + ( pfrld(ji,jj)**rn_betas - pfrld(ji,jj) ) * sprecip(ji,jj)  *             &  ! heat content of precip 
    182                   &      ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus )           & 
    183                   &    + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 ) ) 
    184             ENDIF 
     145            zqld =  tmask(ji,jj,1) * rdt_ice *  & 
     146               &    ( pfrld(ji,jj) * qsr_oce(ji,jj) * frq_m(ji,jj) + pfrld(ji,jj) * qns_oce(ji,jj) + qemp_oce(ji,jj) ) 
    185147 
    186148            ! --- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- ! 
    … …  
    209171            ! Net heat flux on top of ice-ocean [W.m-2] 
    210172            ! ----------------------------------------- 
    211             !     First  step here      : heat flux at the ocean surface + precip 
    212             !     Second step below     : heat flux at the ice   surface (after limthd_dif)  
    213             hfx_in(ji,jj) = hfx_in(ji,jj)                                                                                         &  
    214                ! heat flux above the ocean 
    215                &    +             pfrld(ji,jj)   * ( zqns(ji,jj) + zqsr(ji,jj) )                                                  & 
    216                ! latent heat of precip (note that precip is included in qns but not in qns_ice) 
    217                &    +   ( 1._wp - pfrld(ji,jj) ) * sprecip(ji,jj) * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus )  & 
    218                &    +   ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 ) 
     173            hfx_in(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj)  
    219174 
    220175            ! ----------------------------------------------------------------------------- 
    221             ! Net heat flux that is retroceded to the ocean or taken from the ocean [W.m-2] 
     176            ! Net heat flux on top of the ocean after ice thermo (1st step) [W.m-2] 
    222177            ! ----------------------------------------------------------------------------- 
    223178            !     First  step here              :  non solar + precip - qlead - qturb 
    224179            !     Second step in limthd_dh      :  heat remaining if total melt (zq_rema)  
    225180            !     Third  step in limsbc         :  heat from ice-ocean mass exchange (zf_mass) + solar 
    226             hfx_out(ji,jj) = hfx_out(ji,jj)                                                                                       &  
    227                ! Non solar heat flux received by the ocean 
    228                &    +        pfrld(ji,jj) * qns(ji,jj)                                                                            & 
    229                ! latent heat of precip (note that precip is included in qns but not in qns_ice) 
    230                &    +      ( pfrld(ji,jj)**rn_betas - pfrld(ji,jj) ) * sprecip(ji,jj)       & 
    231                &         * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus )  & 
    232                &    +      ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 )       & 
    233                ! heat flux taken from the ocean where there is open water ice formation 
    234                &    -      qlead(ji,jj) * r1_rdtice                                                                               & 
    235                ! heat flux taken from the ocean during bottom growth/melt (fhld should be 0 while bott growth) 
    236                &    -      at_i(ji,jj) * fhtur(ji,jj)                                                                             & 
    237                &    -      at_i(ji,jj) *  fhld(ji,jj) 
    238  
     181            hfx_out(ji,jj) =   pfrld(ji,jj) * qns_oce(ji,jj) + qemp_oce(ji,jj)  &  ! Non solar heat flux received by the ocean                
     182               &             - qlead(ji,jj) * r1_rdtice                         &  ! heat flux taken from the ocean where there is open water ice formation 
     183               &             - at_i(ji,jj) * fhtur(ji,jj)                       &  ! heat flux taken by turbulence 
     184               &             - at_i(ji,jj) *  fhld(ji,jj)                          ! heat flux taken during bottom growth/melt  
     185                                                                                   !    (fhld should be 0 while bott growth) 
    239186         END DO 
    240187      END DO 
    … …  
    311258            ! --- lateral melting if monocat --- ! 
    312259            !------------------------------------! 
    313             IF ( ( ( nn_monocat == 1 ) .OR. ( nn_monocat == 4 ) ) .AND. ( jpl == 1 ) ) THEN 
     260            IF ( ( nn_monocat == 1 .OR. nn_monocat == 4 ) .AND. jpl == 1 ) THEN 
    314261               CALL lim_thd_lam( 1, nbpb ) 
    315262            END IF 
    … …  
    324271         ENDIF 
    325272         ! 
    326       END DO 
     273      END DO !jl 
    327274 
    328275      !------------------------------------------------------------------------------! 
    … …  
    350297      END DO 
    351298  
    352       !------------------------ 
    353       ! Ice natural aging               
    354       !------------------------ 
    355       oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice /rday 
    356  
    357299      !---------------------------------- 
    358300      ! Change thickness to volume 
    359301      !---------------------------------- 
    360       CALL lim_var_eqv2glo 
     302      v_i(:,:,:)   = ht_i(:,:,:) * a_i(:,:,:) 
     303      v_s(:,:,:)   = ht_s(:,:,:) * a_i(:,:,:) 
     304      smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:) 
     305 
     306      ! update ice age (in case a_i changed, i.e. becomes 0 or lateral melting in monocat) 
     307      DO jl  = 1, jpl 
     308         DO jj = 1, jpj 
     309            DO ji = 1, jpi 
     310               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i_b(ji,jj,jl) - epsi10 ) ) 
     311               oa_i(ji,jj,jl) = rswitch * oa_i(ji,jj,jl) * a_i(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) 
     312            END DO 
     313         END DO 
     314      END DO 
    361315 
    362316      CALL lim_var_zapsmall 
     317 
    363318      !-------------------------------------------- 
    364319      ! Diagnostic thermodynamic growth rates 
    … …  
    399354      ! 
    400355      ! 
    401       CALL wrk_dealloc( jpi, jpj, zqsr, zqns ) 
    402  
    403356      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     357 
    404358      !------------------------------------------------------------------------------| 
    405359      !  6) Transport of ice between thickness categories.                           | 
    406360      !------------------------------------------------------------------------------| 
     361      ! Given thermodynamic growth rates, transport ice between thickness categories. 
    407362      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    408363 
    409       ! Given thermodynamic growth rates, transport ice between thickness categories. 
    410       IF( jpl > 1 )   CALL lim_itd_th_rem( 1, jpl, kt ) 
    411       ! 
    412       CALL lim_var_glo2eqv    ! only for info 
    413       CALL lim_var_agg(1) 
     364      IF( jpl > 1 )      CALL lim_itd_th_rem( 1, jpl, kt ) 
    414365 
    415366      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     367 
    416368      !------------------------------------------------------------------------------| 
    417369      !  7) Add frazil ice growing in leads. 
    418370      !------------------------------------------------------------------------------| 
    419371      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     372 
    420373      CALL lim_thd_lac 
    421       CALL lim_var_glo2eqv    ! only for info 
    422374       
    423       ! conservation test 
    424375      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    425376 
    426       IF(ln_ctl) THEN   ! Control print 
     377      ! Control print 
     378      IF(ln_ctl) THEN 
     379         CALL lim_var_glo2eqv 
     380 
    427381         CALL prt_ctl_info(' ') 
    428382         CALL prt_ctl_info(' - Cell values : ') 
    … …  
    460414   END SUBROUTINE lim_thd  
    461415 
     416  
    462417   SUBROUTINE lim_thd_temp( kideb, kiut ) 
    463418      !!----------------------------------------------------------------------- 
    … …  
    503458      REAL(wp)            ::   zhi_bef            ! ice thickness before thermo 
    504459      REAL(wp)            ::   zdh_mel, zda_mel   ! net melting 
    505       REAL(wp)            ::   zv                 ! ice volume  
     460      REAL(wp)            ::   zvi, zvs           ! ice/snow volumes  
    506461 
    507462      DO ji = kideb, kiut 
    508463         zdh_mel = MIN( 0._wp, dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) ) 
    509          IF( zdh_mel < 0._wp )  THEN 
    510             zv         = a_i_1d(ji) * ht_i_1d(ji) 
     464         IF( zdh_mel < 0._wp .AND. a_i_1d(ji) > 0._wp )  THEN 
     465            zvi          = a_i_1d(ji) * ht_i_1d(ji) 
     466            zvs          = a_i_1d(ji) * ht_s_1d(ji) 
    511467            ! lateral melting = concentration change 
    512468            zhi_bef     = ht_i_1d(ji) - zdh_mel 
    513             zda_mel     =  a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi10 ) ) 
    514             a_i_1d(ji)  = MAX( 0._wp, a_i_1d(ji) + zda_mel )  
    515             ! adjust thickness 
    516             rswitch     = MAX( 0._wp , SIGN( 1._wp , a_i_1d(ji) - epsi20 ) ) 
    517             ht_i_1d(ji) = rswitch * zv / MAX( a_i_1d(ji), epsi20 ) 
     469            rswitch     = MAX( 0._wp , SIGN( 1._wp , zhi_bef - epsi20 ) ) 
     470            zda_mel     = rswitch * a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi20 ) ) 
     471            a_i_1d(ji)  = MAX( epsi20, a_i_1d(ji) + zda_mel )  
     472             ! adjust thickness 
     473            ht_i_1d(ji) = zvi / a_i_1d(ji)             
     474            ht_s_1d(ji) = zvs / a_i_1d(ji)             
    518475            ! retrieve total concentration 
    519476            at_i_1d(ji) = a_i_1d(ji) 
    … …  
    556513         END DO 
    557514          
    558          CALL tab_2d_1d( nbpb, tatm_ice_1d(1:nbpb), tatm_ice(:,:)  , jpi, jpj, npb(1:nbpb) ) 
     515         CALL tab_2d_1d( nbpb, qprec_ice_1d(1:nbpb), qprec_ice(:,:) , jpi, jpj, npb(1:nbpb) ) 
    559516         CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb), qsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 
    560517         CALL tab_2d_1d( nbpb, fr1_i0_1d  (1:nbpb), fr1_i0          , jpi, jpj, npb(1:nbpb) ) 
    … …  
    562519         CALL tab_2d_1d( nbpb, qns_ice_1d (1:nbpb), qns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 
    563520         CALL tab_2d_1d( nbpb, ftr_ice_1d (1:nbpb), ftr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 
    564          IF( .NOT. lk_cpl ) THEN 
    565             CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb), qla_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 
    566             CALL tab_2d_1d( nbpb, dqla_ice_1d(1:nbpb), dqla_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 
    567          ENDIF 
     521         CALL tab_2d_1d( nbpb, evap_ice_1d (1:nbpb), evap_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 
    568522         CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 
    569523         CALL tab_2d_1d( nbpb, t_bo_1d     (1:nbpb), t_bo            , jpi, jpj, npb(1:nbpb) ) 
    … …  
    656610         CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj) 
    657611         CALL tab_1d_2d( nbpb, ftr_ice(:,:,jl), npb, ftr_ice_1d(1:nbpb) , jpi, jpj ) 
    658                    
     612         !          
    659613      END SELECT 
    660614 
    … …  
    676630      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    677631      NAMELIST/namicethd/ rn_hnewice, ln_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb,                       & 
    678          &                rn_himin, parsub, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, & 
     632         &                rn_himin, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, & 
    679633         &                nn_monocat, ln_it_qnsice 
    680634      !!------------------------------------------------------------------- 
    … …  
    700654      ENDIF 
    701655 
    702       IF( lk_cpl .AND. parsub /= 0.0 )   CALL ctl_stop( 'In coupled mode, use parsub = 0. or send dqla' ) 
    703656      ! 
    704657      IF(lwp) THEN                          ! control print 
    … …  
    712665         WRITE(numout,*)'      minimum ice thickness                                   rn_himin     = ', rn_himin  
    713666         WRITE(numout,*)'      numerical carac. of the scheme for diffusion in ice ' 
    714          WRITE(numout,*)'      switch for snow sublimation  (=1) or not (=0)           parsub       = ', parsub   
    715667         WRITE(numout,*)'      coefficient for ice-lead partition of snowfall          rn_betas     = ', rn_betas 
    716668         WRITE(numout,*)'      extinction radiation parameter in sea ice               rn_kappa_i   = ', rn_kappa_i 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

    r5134 r5621  
    2929   PRIVATE 
    3030 
    31    PUBLIC   lim_thd_dh   ! called by lim_thd 
     31   PUBLIC   lim_thd_dh      ! called by lim_thd 
     32   PUBLIC   lim_thd_snwblow ! called in sbcblk/sbcclio/sbccpl and here 
     33 
     34   INTERFACE lim_thd_snwblow 
     35      MODULE PROCEDURE lim_thd_snwblow_1d, lim_thd_snwblow_2d 
     36   END INTERFACE 
    3237 
    3338   !!---------------------------------------------------------------------- 
    … …  
    7176      REAL(wp) ::   zfdum        
    7277      REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment 
    73       REAL(wp) ::   zcoeff       ! dummy argument for snowfall partitioning over ice and leads 
    74       REAL(wp) ::   zs_snic  ! snow-ice salinity 
     78      REAL(wp) ::   zs_snic      ! snow-ice salinity 
    7579      REAL(wp) ::   zswi1        ! switch for computation of bottom salinity 
    7680      REAL(wp) ::   zswi12       ! switch for computation of bottom salinity 
    … …  
    8690      REAL(wp) ::   zsstK        ! SST in Kelvin 
    8791 
    88       REAL(wp), POINTER, DIMENSION(:) ::   zh_s        ! snow layer thickness 
    8992      REAL(wp), POINTER, DIMENSION(:) ::   zqprec      ! energy of fallen snow                       (J.m-3) 
    9093      REAL(wp), POINTER, DIMENSION(:) ::   zq_su       ! heat for surface ablation                   (J.m-2) 
    … …  
    9295      REAL(wp), POINTER, DIMENSION(:) ::   zq_rema     ! remaining heat at the end of the routine    (J.m-2) 
    9396      REAL(wp), POINTER, DIMENSION(:) ::   zf_tt       ! Heat budget to determine melting or freezing(W.m-2) 
    94       INTEGER , POINTER, DIMENSION(:) ::   icount      ! number of layers vanished by melting  
    9597 
    9698      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_mel   ! snow melt  
    … …  
    100102      REAL(wp), POINTER, DIMENSION(:,:) ::   zdeltah 
    101103      REAL(wp), POINTER, DIMENSION(:,:) ::   zh_i      ! ice layer thickness 
     104      INTEGER , POINTER, DIMENSION(:,:) ::   icount    ! number of layers vanished by melting  
    102105 
    103106      REAL(wp), POINTER, DIMENSION(:) ::   zqh_i       ! total ice heat content  (J.m-2) 
    104107      REAL(wp), POINTER, DIMENSION(:) ::   zqh_s       ! total snow heat content (J.m-2) 
    105108      REAL(wp), POINTER, DIMENSION(:) ::   zq_s        ! total snow enthalpy     (J.m-3) 
     109      REAL(wp), POINTER, DIMENSION(:) ::   zsnw        ! distribution of snow after wind blowing 
    106110 
    107111      REAL(wp) :: zswitch_sal 
    … …  
    118122      END SELECT 
    119123 
    120       CALL wrk_alloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) 
     124      CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw ) 
    121125      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 
    122       CALL wrk_alloc( jpij, nlay_i+1, zdeltah, zh_i ) 
    123       CALL wrk_alloc( jpij, icount ) 
    124        
     126      CALL wrk_alloc( jpij, nlay_i, zdeltah, zh_i ) 
     127      CALL wrk_alloc( jpij, nlay_i, icount ) 
     128        
    125129      dh_i_surf  (:) = 0._wp ; dh_i_bott  (:) = 0._wp ; dh_snowice(:) = 0._wp 
    126130      dsm_i_se_1d(:) = 0._wp ; dsm_i_si_1d(:) = 0._wp    
    127   
    128       zqprec (:) = 0._wp ; zq_su  (:) = 0._wp ; zq_bo  (:) = 0._wp ; zf_tt  (:) = 0._wp 
    129       zq_rema(:) = 0._wp 
    130  
    131       zh_s     (:) = 0._wp        
    132       zdh_s_pre(:) = 0._wp 
    133       zdh_s_mel(:) = 0._wp 
    134       zdh_s_sub(:) = 0._wp 
    135       zqh_s    (:) = 0._wp       
    136       zqh_i    (:) = 0._wp    
    137  
    138       zh_i      (:,:) = 0._wp        
    139       zdeltah   (:,:) = 0._wp        
    140       icount    (:)   = 0 
     131 
     132      zqprec   (:) = 0._wp ; zq_su    (:) = 0._wp ; zq_bo    (:) = 0._wp ; zf_tt(:) = 0._wp 
     133      zq_rema  (:) = 0._wp ; zsnw     (:) = 0._wp 
     134      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      
     136 
     137      zdeltah(:,:) = 0._wp ; zh_i(:,:) = 0._wp        
     138      icount (:,:) = 0 
     139 
     140 
     141      ! Initialize enthalpy at nlay_i+1 
     142      DO ji = kideb, kiut 
     143         q_i_1d(ji,nlay_i+1) = 0._wp 
     144      END DO 
    141145 
    142146      ! initialize layer thicknesses and enthalpies 
    … …  
    155159      ! 
    156160      DO ji = kideb, kiut 
    157          rswitch       = 1._wp - MAX(  0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) 
    158          ztmelts       = rswitch * rt0 + ( 1._wp - rswitch ) * rt0 
    159  
    160161         zfdum      = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji)  
    161162         zf_tt(ji)  = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji)  
    162163 
    163          zq_su (ji) = MAX( 0._wp, zfdum     * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - ztmelts ) ) 
     164         zq_su (ji) = MAX( 0._wp, zfdum     * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) ) 
    164165         zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice ) 
    165166      END DO 
    … …  
    187188      !------------------------------------------------------------! 
    188189      ! 
    189       DO ji = kideb, kiut      
    190          zh_s(ji) = ht_s_1d(ji) * r1_nlay_s 
    191       END DO 
    192       ! 
    193190      DO jk = 1, nlay_s 
    194191         DO ji = kideb, kiut 
    195             zqh_s(ji) =  zqh_s(ji) + q_s_1d(ji,jk) * zh_s(ji) 
     192            zqh_s(ji) =  zqh_s(ji) + q_s_1d(ji,jk) * ht_s_1d(ji) * r1_nlay_s 
    196193         END DO 
    197194      END DO 
    … …  
    222219      ! Martin Vancoppenolle, December 2006 
    223220 
     221      CALL lim_thd_snwblow( 1. - at_i_1d(kideb:kiut), zsnw(kideb:kiut) ) ! snow distribution over ice after wind blowing 
     222 
     223      zdeltah(:,:) = 0._wp 
    224224      DO ji = kideb, kiut 
    225225         !----------- 
    … …  
    227227         !----------- 
    228228         ! thickness change 
    229          zcoeff = ( 1._wp - ( 1._wp - at_i_1d(ji) )**rn_betas ) / at_i_1d(ji)  
    230          zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice * r1_rhosn 
    231          ! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K) 
    232          zqprec   (ji) = rhosn * ( cpic * ( rt0 - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus )    
     229         zdh_s_pre(ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhosn / at_i_1d(ji) 
     230         ! enthalpy of the precip (>0, J.m-3) 
     231         zqprec   (ji) = - qprec_ice_1d(ji)    
    233232         IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp 
    234233         ! heat flux from snow precip (>0, W.m-2) 
    … …  
    236235         ! mass flux, <0 
    237236         wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice 
    238          ! update thickness 
    239          ht_s_1d    (ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) ) 
    240237 
    241238         !--------------------- 
    … …  
    243240         !--------------------- 
    244241         ! thickness change 
    245          IF( zdh_s_pre(ji) > 0._wp ) THEN 
    246242         rswitch        = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) ) 
    247          zdh_s_mel (ji) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 ) 
    248          zdh_s_mel (ji) = MAX( - zdh_s_pre(ji), zdh_s_mel(ji) ) ! bound melting  
     243         zdeltah (ji,1) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 ) 
     244         zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting  
    249245         ! heat used to melt snow (W.m-2, >0) 
    250          hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdh_s_mel(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice 
     246         hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_rdtice 
    251247         ! snow melting only = water into the ocean (then without snow precip), >0 
    252          wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_mel(ji) * r1_rdtice 
    253           
    254          ! updates available heat + thickness 
    255          zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdh_s_mel(ji) * zqprec(ji) )       
    256          ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_mel(ji) ) 
    257          zh_s  (ji) = ht_s_1d(ji) * r1_nlay_s 
    258  
    259          ENDIF 
    260       END DO 
    261  
    262       ! If heat still available, then melt more snow 
    263       zdeltah(:,:) = 0._wp ! important 
     248         wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice     
     249         ! updates available heat + precipitations after melting 
     250         zq_su     (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) )       
     251         zdh_s_pre (ji) = zdh_s_pre(ji) + zdeltah(ji,1) 
     252 
     253         ! update thickness 
     254         ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) ) 
     255      END DO 
     256 
     257      ! If heat still available (zq_su > 0), then melt more snow 
     258      zdeltah(:,:) = 0._wp 
    264259      DO jk = 1, nlay_s 
    265260         DO ji = kideb, kiut 
    … …  
    268263            rswitch          = rswitch * ( MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,jk) - epsi20 ) ) )  
    269264            zdeltah  (ji,jk) = - rswitch * zq_su(ji) / MAX( q_s_1d(ji,jk), epsi20 ) 
    270             zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) ! bound melting 
     265            zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - ht_s_1d(ji) ) ! bound melting 
    271266            zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk)     
    272267            ! heat used to melt snow(W.m-2, >0) 
    … …  
    274269            ! snow melting only = water into the ocean (then without snow precip) 
    275270            wfx_snw_1d(ji)   = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
    276  
    277271            ! updates available heat + thickness 
    278272            zq_su (ji)  = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) ) 
    279273            ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) ) 
    280  
    281274         END DO 
    282275      END DO 
    … …  
    286279      !---------------------- 
    287280      ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates 
    288       ! clem comment: not counted in mass exchange in limsbc since this is an exchange with atm. (not ocean) 
     281      ! clem comment: not counted in mass/heat exchange in limsbc since this is an exchange with atm. (not ocean) 
    289282      ! clem comment: ice should also sublimate 
    290       IF( lk_cpl ) THEN 
    291          ! coupled mode: sublimation already included in emp_ice (to do in limsbc_ice) 
    292          zdh_s_sub(:)      =  0._wp  
    293       ELSE 
    294          ! forced  mode: snow thickness change due to sublimation 
    295          DO ji = kideb, kiut 
    296             zdh_s_sub(ji)  =  MAX( - ht_s_1d(ji) , - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice ) 
    297             ! Heat flux by sublimation [W.m-2], < 0 
    298             !      sublimate first snow that had fallen, then pre-existing snow 
    299             zcoeff         =      ( MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) )   * zqprec(ji) +   & 
    300                &  ( zdh_s_sub(ji) - MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) ) * q_s_1d(ji,1) )  & 
    301                &  * a_i_1d(ji) * r1_rdtice 
    302             hfx_sub_1d(ji) = hfx_sub_1d(ji) + zcoeff 
    303             ! Mass flux by sublimation 
    304             wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice 
    305             ! new snow thickness 
    306             ht_s_1d(ji)     =  MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) ) 
    307          END DO 
    308       ENDIF 
    309  
     283      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 
     290         zdeltah(ji,1)  = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 
     291         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)  & 
     292            &                              ) * a_i_1d(ji) * r1_rdtice 
     293         ! Mass flux by sublimation 
     294         wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice 
     295         ! new snow thickness 
     296         ht_s_1d(ji)    =  MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) ) 
     297         ! update precipitations after sublimation and correct sublimation 
     298         zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1) 
     299         zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1) 
     300      END DO 
     301       
    310302      ! --- Update snow diags --- ! 
    311303      DO ji = kideb, kiut 
    312          dh_s_tot(ji)   = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 
    313          zh_s(ji)       = ht_s_1d(ji) * r1_nlay_s 
    314       END DO ! ji 
     304         dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 
     305      END DO 
    315306 
    316307      !------------------------------------------- 
    … …  
    323314            rswitch       = MAX(  0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 )  ) 
    324315            q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) *                          & 
    325               &            ( (   MAX( 0._wp, dh_s_tot(ji) )               ) * zqprec(ji) +  & 
    326               &              ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_1d(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) ) 
     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 ) ) 
    327318            zq_s(ji)     =  zq_s(ji) + q_s_1d(ji,jk) 
    328319         END DO 
    … …  
    334325      zdeltah(:,:) = 0._wp ! important 
    335326      DO jk = 1, nlay_i 
    336          DO ji = kideb, kiut  
    337             zEi            = - q_i_1d(ji,jk) * r1_rhoic             ! Specific enthalpy of layer k [J/kg, <0] 
    338  
    339             ztmelts        = - tmut * s_i_1d(ji,jk) + rt0           ! Melting point of layer k [K] 
    340  
    341             zEw            =    rcp * ( ztmelts - rt0 )            ! Specific enthalpy of resulting meltwater [J/kg, <0] 
    342  
    343             zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0 
    344  
    345             zfmdt          = - zq_su(ji) / zdE                     ! Mass flux to the ocean [kg/m2, >0] 
    346  
    347             zdeltah(ji,jk) = - zfmdt * r1_rhoic                    ! Melt of layer jk [m, <0] 
    348  
    349             zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0] 
    350  
    351             zq_su(ji)      = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat 
    352  
    353             dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt 
    354  
    355             zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0] 
    356  
    357             zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0] 
    358  
    359             ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 
    360             sfx_sum_1d(ji)   = sfx_sum_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice 
    361  
    362             ! Contribution to heat flux [W.m-2], < 0 
    363             hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
    364  
    365             ! Total heat flux used in this process [W.m-2], > 0   
    366             hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 
    367  
    368             ! Contribution to mass flux 
    369             wfx_sum_1d(ji) =  wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
    370             
     327         DO ji = kideb, kiut 
     328            ztmelts           = - tmut * s_i_1d(ji,jk) + rt0          ! Melting point of layer k [K] 
     329             
     330            IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 
     331 
     332               zEi            = - q_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]        
     333               zdE            = 0._wp                                 ! Specific enthalpy difference   (J/kg, <0) 
     334                                                                      ! set up at 0 since no energy is needed to melt water...(it is already melted) 
     335               zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )          ! internal melting occurs when the internal temperature is above freezing      
     336                                                                      ! this should normally not happen, but sometimes, heat diffusion leads to this 
     337               zfmdt          = - zdeltah(ji,jk) * rhoic              ! Mass flux x time step > 0 
     338                          
     339               dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt 
     340                
     341               zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0] 
     342 
     343               ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)  
     344               hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice 
     345                
     346               ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 
     347               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 
     348                
     349               ! Contribution to mass flux 
     350               wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
     351 
     352            ELSE                                !!! Surface melting 
     353                
     354               zEi            = - q_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0] 
     355               zEw            =    rcp * ( ztmelts - rt0 )            ! Specific enthalpy of resulting meltwater [J/kg, <0] 
     356               zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0 
     357                
     358               zfmdt          = - zq_su(ji) / zdE                     ! Mass flux to the ocean [kg/m2, >0] 
     359                
     360               zdeltah(ji,jk) = - zfmdt * r1_rhoic                    ! Melt of layer jk [m, <0] 
     361                
     362               zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0] 
     363                
     364               zq_su(ji)      = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat 
     365                
     366               dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt 
     367                
     368               zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0] 
     369                
     370               zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0] 
     371                
     372               ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 
     373               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 
     374                
     375               ! Contribution to heat flux [W.m-2], < 0 
     376               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
     377                
     378               ! Total heat flux used in this process [W.m-2], > 0   
     379               hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 
     380                
     381               ! Contribution to mass flux 
     382               wfx_sum_1d(ji) =  wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
     383                
     384            END IF 
    371385            ! record which layers have disappeared (for bottom melting)  
    372386            !    => icount=0 : no layer has vanished 
    373387            !    => icount=5 : 5 layers have vanished 
    374             rswitch     = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) )  
    375             icount(ji)  = icount(ji) + NINT( rswitch ) 
    376             zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) ) 
     388            rswitch       = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) )  
     389            icount(ji,jk) = NINT( rswitch ) 
     390            zh_i(ji,jk)   = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) ) 
    377391 
    378392            ! update heat content (J.m-2) and layer thickness 
    … …  
    405419      ! -> need for an iterative procedure, which converges quickly 
    406420 
    407       IF ( nn_icesal == 2 ) THEN 
    408          num_iter_max = 5 
    409       ELSE 
    410          num_iter_max = 1 
    411       ENDIF 
    412  
    413       ! Just to be sure that enthalpy at nlay_i+1 is null 
    414       DO ji = kideb, kiut 
    415          q_i_1d(ji,nlay_i+1) = 0._wp 
    416       END DO 
     421      num_iter_max = 1 
     422      IF( nn_icesal == 2 ) num_iter_max = 5 
    417423 
    418424      ! Iterative procedure 
    … …  
    483489             
    484490            ! Contribution to salt flux, <0 
    485             sfx_bog_1d(ji) = sfx_bog_1d(ji) + s_i_new(ji) * a_i_1d(ji) * zfmdt * r1_rdtice 
     491            sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * s_i_new(ji) * r1_rdtice 
    486492 
    487493            ! Contribution to mass flux, <0 
    … …  
    500506      DO jk = nlay_i, 1, -1 
    501507         DO ji = kideb, kiut 
    502             IF(  zf_tt(ji)  >=  0._wp  .AND. jk > icount(ji) ) THEN   ! do not calculate where layer has already disappeared by surface melting  
     508            IF(  zf_tt(ji)  >  0._wp  .AND. jk > icount(ji,jk) ) THEN   ! do not calculate where layer has already disappeared by surface melting  
    503509 
    504510               ztmelts = - tmut * s_i_1d(ji,jk) + rt0  ! Melting point of layer jk (K) 
    … …  
    507513 
    508514                  zEi               = - q_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0) 
    509  
    510                   !!zEw               = rcp * ( t_i_1d(ji,jk) - rt0 )  ! Specific enthalpy of meltwater at T = t_i_1d (J/kg, <0) 
    511  
    512515                  zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0) 
    513516                                                                    ! set up at 0 since no energy is needed to melt water...(it is already melted) 
    514  
    515                   zdeltah   (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing      
    516                                                                    ! this should normally not happen, but sometimes, heat diffusion leads to this 
     517                  zdeltah   (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )  ! internal melting occurs when the internal temperature is above freezing      
     518                                                                    ! this should normally not happen, but sometimes, heat diffusion leads to this 
    517519 
    518520                  dh_i_bott (ji)    = dh_i_bott(ji) + zdeltah(ji,jk) 
    519521 
    520                   zfmdt             = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step > 0 
     522                  zfmdt             = - zdeltah(ji,jk) * rhoic      ! Mass flux x time step > 0 
    521523 
    522524                  ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)  
    … …  
    524526 
    525527                  ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 
    526                   sfx_res_1d(ji) = sfx_res_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice 
     528                  sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 
    527529                                     
    528530                  ! Contribution to mass flux 
    … …  
    535537               ELSE                               !!! Basal melting 
    536538 
    537                   zEi               = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 
    538  
    539                   zEw               = rcp * ( ztmelts - rt0 )    ! Specific enthalpy of meltwater (J/kg, <0) 
    540  
    541                   zdE               = zEi - zEw                  ! Specific enthalpy difference   (J/kg, <0) 
    542  
    543                   zfmdt             = - zq_bo(ji) / zdE          ! Mass flux x time step (kg/m2, >0) 
    544  
    545                   zdeltah(ji,jk)    = - zfmdt * r1_rhoic         ! Gross thickness change 
    546  
    547                   zdeltah(ji,jk)    = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change 
     539                  zEi             = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 
     540                  zEw             = rcp * ( ztmelts - rt0 )    ! Specific enthalpy of meltwater (J/kg, <0) 
     541                  zdE             = zEi - zEw                  ! Specific enthalpy difference   (J/kg, <0) 
     542 
     543                  zfmdt           = - zq_bo(ji) / zdE          ! Mass flux x time step (kg/m2, >0) 
     544 
     545                  zdeltah(ji,jk)  = - zfmdt * r1_rhoic         ! Gross thickness change 
     546 
     547                  zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change 
    548548                   
    549                   zq_bo(ji)         = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors 
    550  
    551                   dh_i_bott(ji)     = dh_i_bott(ji) + zdeltah(ji,jk)    ! Update basal melt 
    552  
    553                   zfmdt             = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step > 0 
    554  
    555                   zQm               = zfmdt * zEw         ! Heat exchanged with ocean 
     549                  zq_bo(ji)       = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors 
     550 
     551                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)    ! Update basal melt 
     552 
     553                  zfmdt           = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step > 0 
     554 
     555                  zQm             = zfmdt * zEw         ! Heat exchanged with ocean 
    556556 
    557557                  ! Contribution to heat flux to the ocean [W.m-2], <0   
    558                   hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
     558                  hfx_thd_1d(ji)  = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
    559559 
    560560                  ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 
    561                   sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice 
     561                  sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoic *  a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 
    562562                   
    563563                  ! Total heat flux used in this process [W.m-2], >0   
    564                   hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 
     564                  hfx_bom_1d(ji)  = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 
    565565                   
    566566                  ! Contribution to mass flux 
    567                   wfx_bom_1d(ji) =  wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
     567                  wfx_bom_1d(ji)  =  wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
    568568 
    569569                  ! update heat content (J.m-2) and layer thickness 
    … …  
    595595         zdeltah  (ji,1) = - rswitch * zq_rema(ji) / MAX( q_s_1d(ji,1), epsi20 ) 
    596596         zdeltah  (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting 
    597          zdh_s_mel(ji)   = zdh_s_mel(ji) + zdeltah(ji,1)     
    598597         dh_s_tot (ji)   = dh_s_tot(ji)  + zdeltah(ji,1) 
    599598         ht_s_1d   (ji)  = ht_s_1d(ji)   + zdeltah(ji,1) 
    … …  
    622621         dh_snowice(ji) = MAX(  0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) / ( rhosn+rau0-rhoic )  ) 
    623622 
    624          ht_i_1d(ji)     = ht_i_1d(ji) + dh_snowice(ji) 
    625          ht_s_1d(ji)     = ht_s_1d(ji) - dh_snowice(ji) 
     623         ht_i_1d(ji)    = ht_i_1d(ji) + dh_snowice(ji) 
     624         ht_s_1d(ji)    = ht_s_1d(ji) - dh_snowice(ji) 
    626625 
    627626         ! Salinity of snow ice 
    … …  
    669668      ! Update temperature, energy 
    670669      !------------------------------------------- 
    671       !clem bug: we should take snow into account here 
    672670      DO ji = kideb, kiut 
    673671         rswitch     =  1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) )  
    … …  
    688686      WHERE( ht_i_1d == 0._wp ) a_i_1d = 0._wp 
    689687       
    690       CALL wrk_dealloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) 
     688      CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw ) 
    691689      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 
    692       CALL wrk_dealloc( jpij, nlay_i+1, zdeltah, zh_i ) 
    693       CALL wrk_dealloc( jpij, icount ) 
     690      CALL wrk_dealloc( jpij, nlay_i, zdeltah, zh_i ) 
     691      CALL wrk_dealloc( jpij, nlay_i, icount ) 
    694692      ! 
    695693      ! 
    696694   END SUBROUTINE lim_thd_dh 
     695 
     696 
     697   !!-------------------------------------------------------------------------- 
     698   !! INTERFACE lim_thd_snwblow 
     699   !! ** Purpose :   Compute distribution of precip over the ice 
     700   !!-------------------------------------------------------------------------- 
     701   SUBROUTINE lim_thd_snwblow_2d( pin, pout ) 
     702      REAL(wp), DIMENSION(:,:), INTENT(in   ) :: pin   ! previous fraction lead ( pfrld or (1. - a_i_b) ) 
     703      REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout 
     704      pout = ( 1._wp - ( pin )**rn_betas ) 
     705   END SUBROUTINE lim_thd_snwblow_2d 
     706 
     707   SUBROUTINE lim_thd_snwblow_1d( pin, pout ) 
     708      REAL(wp), DIMENSION(:), INTENT(in   ) :: pin 
     709      REAL(wp), DIMENSION(:), INTENT(inout) :: pout 
     710      pout = ( 1._wp - ( pin )**rn_betas ) 
     711   END SUBROUTINE lim_thd_snwblow_1d 
     712 
    697713    
    698714#else 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90

    r5146 r5621  
    2424   USE wrk_nemo       ! work arrays 
    2525   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    26    USE sbc_oce, ONLY : lk_cpl 
    2726 
    2827   IMPLICIT NONE 
    … …  
    170169      CALL wrk_alloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw ) 
    171170      CALL wrk_alloc( jpij, zf, dzf, zqns_ice_b, zerrit, zdifcase, zftrice, zihic, zghe ) 
    172       CALL wrk_alloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0) 
    173       CALL wrk_alloc( jpij, nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0) 
    174       CALL wrk_alloc( jpij, nlay_i+3, zindterm, zindtbis, zdiagbis  ) 
    175       CALL wrk_alloc( jpij, nlay_i+3, 3, ztrid ) 
     171      CALL wrk_alloc( jpij,nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0 ) 
     172      CALL wrk_alloc( jpij,nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0 ) 
     173      CALL wrk_alloc( jpij,nlay_i+3, zindterm, zindtbis, zdiagbis  ) 
     174      CALL wrk_alloc( jpij,nlay_i+3,3, ztrid ) 
    176175 
    177176      CALL wrk_alloc( jpij, zdq, zq_ini, zhfx_err ) 
    … …  
    283282      END DO 
    284283 
    285       ! 
    286284      !------------------------------------------------------------------------------| 
    287285      !  3) Iterative procedure begins                                               | 
    … …  
    679677         END DO 
    680678 
    681          DO numeq = nlay_i + nlay_s + 1, nlay_s + 2, -1 
     679         DO numeq = nlay_i + nlay_s, nlay_s + 2, -1 
    682680            DO ji = kideb , kiut 
    683681               jk    =  numeq - nlay_s - 1 
    … …  
    746744      !-------------------------------------------------------------------------! 
    747745      DO ji = kideb, kiut 
    748          ! forced mode only : update of latent heat fluxes (sublimation) (always >=0, upward flux)  
    749          IF( .NOT. lk_cpl) qla_ice_1d (ji) = MAX( 0._wp, qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) ) 
    750746         !                                ! surface ice conduction flux 
    751747         isnow(ji)       = 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) ) 
    … …  
    760756 
    761757      ! --- diagnose the change in non-solar flux due to surface temperature change --- ! 
    762       IF ( ln_it_qnsice ) hfx_err_dif_1d(:) = hfx_err_dif_1d(:) - ( qns_ice_1d(:)  - zqns_ice_b(:) ) * a_i_1d(:)  
     758      IF ( ln_it_qnsice ) THEN 
     759         DO ji = kideb, kiut 
     760            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji)  - zqns_ice_b(ji) ) * a_i_1d(ji)  
     761         END DO 
     762      END IF 
    763763 
    764764      ! --- diag conservation imbalance on heat diffusion - PART 2 --- ! 
    … …  
    772772         ENDIF 
    773773         hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_1d(ji) 
     774 
     775         ! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation) 
     776         hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err(ji) * a_i_1d(ji) 
    774777      END DO  
    775  
    776       hfx_err_dif_1d(:) = hfx_err_dif_1d(:) - zhfx_err(:) * a_i_1d(:) 
    777778 
    778779      !----------------------------------------- 
    … …  
    787788               &             ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 
    788789         ENDIF 
    789       END DO 
    790  
    791       ! --- compute diagnostic net heat flux at the surface of the snow-ice system (W.m-2) 
    792       DO ji = kideb, kiut 
    793          ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 
    794          hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_1d(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) ) 
    795       END DO 
    796     
     790         ! correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 
     791         hfx_dif_1d(ji) = hfx_dif_1d(ji) - zhfx_err(ji) * a_i_1d(ji) 
     792      END DO    
    797793      ! 
    798794      CALL wrk_dealloc( jpij, numeqmin, numeqmax ) 
    799795      CALL wrk_dealloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw ) 
    800       CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zghe ) 
    801       CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i,   & 
    802          &              ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 
    803       CALL wrk_dealloc( jpij, nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 ) 
    804       CALL wrk_dealloc( jpij, nlay_i+3, zindterm, zindtbis, zdiagbis ) 
    805       CALL wrk_dealloc( jpij, nlay_i+3, 3, ztrid ) 
     796      CALL wrk_dealloc( jpij, zf, dzf, zqns_ice_b, zerrit, zdifcase, zftrice, zihic, zghe ) 
     797      CALL wrk_dealloc( jpij,nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 
     798      CALL wrk_dealloc( jpij,nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 ) 
     799      CALL wrk_dealloc( jpij,nlay_i+3, zindterm, zindtbis, zdiagbis ) 
     800      CALL wrk_dealloc( jpij,nlay_i+3,3, ztrid ) 
    806801      CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx_err ) 
    807802 
    … …  
    824819      DO jk = 1, nlay_i             ! Sea ice energy of melting 
    825820         DO ji = kideb, kiut 
    826             ztmelts      = - tmut  * s_i_1d(ji,jk) + rt0  
    827             rswitch      = MAX( 0._wp , SIGN( 1._wp , -(t_i_1d(ji,jk) - rt0) - epsi20 ) ) 
    828             q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) )                                                & 
    829                &                   + lfus * ( 1.0 - rswitch * ( ztmelts-rt0 ) / MIN( t_i_1d(ji,jk) - rt0, -epsi20 ) )   & 
    830                &                   - rcp  *                   ( ztmelts-rt0 )  )  
     821            ztmelts      = - tmut  * s_i_1d(ji,jk) + rt0 
     822            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts ) ! Force t_i_1d to be lower than melting point 
     823                                                          !   (sometimes dif scheme produces abnormally high temperatures)    
     824            q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) )                           & 
     825               &                    + lfus * ( 1.0 - ( ztmelts-rt0 ) / ( t_i_1d(ji,jk) - rt0 ) )   & 
     826               &                    - rcp  *         ( ztmelts-rt0 )  )  
    831827         END DO 
    832828      END DO 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90

    r5134 r5621  
    3131   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3232   USE limthd_ent 
     33   USE limvar 
    3334 
    3435   IMPLICIT NONE 
    … …  
    105106      REAL(wp), POINTER, DIMENSION(:,:) ::   za_i_1d   ! 1-D version of a_i 
    106107      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_i_1d   ! 1-D version of v_i 
    107       REAL(wp), POINTER, DIMENSION(:,:) ::   zoa_i_1d  ! 1-D version of oa_i 
    108108      REAL(wp), POINTER, DIMENSION(:,:) ::   zsmv_i_1d ! 1-D version of smv_i 
    109109 
    … …  
    118118      CALL wrk_alloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 
    119119      CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d ) 
    120       CALL wrk_alloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 
    121       CALL wrk_alloc( jpij,nlay_i+1,jpl, ze_i_1d ) 
     120      CALL wrk_alloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zsmv_i_1d ) 
     121      CALL wrk_alloc( jpij,nlay_i,jpl, ze_i_1d ) 
    122122      CALL wrk_alloc( jpi,jpj, zvrel ) 
    123123 
     124      CALL lim_var_agg(1) 
     125      CALL lim_var_glo2eqv 
    124126      !------------------------------------------------------------------------------| 
    125127      ! 2) Convert units for ice internal energy 
    … …  
    289291            CALL tab_2d_1d( nbpac, za_i_1d  (1:nbpac,jl), a_i  (:,:,jl), jpi, jpj, npac(1:nbpac) ) 
    290292            CALL tab_2d_1d( nbpac, zv_i_1d  (1:nbpac,jl), v_i  (:,:,jl), jpi, jpj, npac(1:nbpac) ) 
    291             CALL tab_2d_1d( nbpac, zoa_i_1d (1:nbpac,jl), oa_i (:,:,jl), jpi, jpj, npac(1:nbpac) ) 
    292293            CALL tab_2d_1d( nbpac, zsmv_i_1d(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) ) 
    293294            DO jk = 1, nlay_i 
    294295               CALL tab_2d_1d( nbpac, ze_i_1d(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) ) 
    295             END DO ! jk 
    296          END DO ! jl 
     296            END DO 
     297         END DO 
    297298 
    298299         CALL tab_2d_1d( nbpac, qlead_1d  (1:nbpac)     , qlead  , jpi, jpj, npac(1:nbpac) ) 
    … …  
    355356         DO ji = 1, nbpac 
    356357            zo_newice(ji) = 0._wp 
    357          END DO ! ji 
     358         END DO 
    358359 
    359360         !------------------- 
    … …  
    477478         ENDDO 
    478479 
    479          !------------ 
    480          ! Update age  
    481          !------------ 
    482          DO jl = 1, jpl 
    483             DO ji = 1, nbpac 
    484                rswitch          = MAX( 0._wp , SIGN( 1._wp , za_i_1d(ji,jl) - epsi20 ) )  ! 0 if no ice and 1 if yes 
    485                zoa_i_1d(ji,jl)  = za_b(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * rswitch    
    486             END DO  
    487          END DO    
    488  
    489480         !----------------- 
    490481         ! Update salinity 
    … …  
    503494            CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_1d (1:nbpac,jl), jpi, jpj ) 
    504495            CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_1d (1:nbpac,jl), jpi, jpj ) 
    505             CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_1d(1:nbpac,jl), jpi, jpj ) 
    506496            CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_1d(1:nbpac,jl) , jpi, jpj ) 
    507497            DO jk = 1, nlay_i 
    … …  
    535525      CALL wrk_dealloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 
    536526      CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d ) 
    537       CALL wrk_dealloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 
    538       CALL wrk_dealloc( jpij,nlay_i+1,jpl, ze_i_1d ) 
     527      CALL wrk_dealloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zsmv_i_1d ) 
     528      CALL wrk_dealloc( jpij,nlay_i,jpl, ze_i_1d ) 
    539529      CALL wrk_dealloc( jpi,jpj, zvrel ) 
    540530      ! 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90

    r5138 r5621  
    8080      IF( nn_timing == 1 )  CALL timing_start('limtrp') 
    8181 
    82       CALL wrk_alloc( jpi,jpj,           zsm, zatold, zeiold, zesold ) 
    83       CALL wrk_alloc( jpi,jpj,jpl,       z0ice, z0snw, z0ai, z0es , z0smi , z0oi ) 
    84       CALL wrk_alloc( jpi,jpj,1,         z0opw ) 
    85       CALL wrk_alloc( jpi,jpj,nlay_i+1,jpl, z0ei ) 
    86       CALL wrk_alloc( jpi,jpj,jpl,       zhimax, zviold, zvsold, zsmvold ) 
     82      CALL wrk_alloc( jpi,jpj,            zsm, zatold, zeiold, zesold ) 
     83      CALL wrk_alloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi ) 
     84      CALL wrk_alloc( jpi,jpj,1,          z0opw ) 
     85      CALL wrk_alloc( jpi,jpj,nlay_i,jpl, z0ei ) 
     86      CALL wrk_alloc( jpi,jpj,jpl,        zhimax, zviold, zvsold, zsmvold ) 
    8787 
    8888      IF( numit == nstart .AND. lwp ) THEN 
    … …  
    112112 
    113113         !--- Thickness correction init. ------------------------------- 
    114          CALL lim_var_glo2eqv 
    115114         zatold(:,:) = SUM( a_i(:,:,:), dim=3 ) 
     115         DO jl = 1, jpl 
     116            DO jj = 1, jpj 
     117               DO ji = 1, jpi 
     118                  rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) 
     119                  ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 
     120                  ht_s  (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 
     121               END DO 
     122            END DO 
     123         END DO 
    116124         !--------------------------------------------------------------------- 
    117          ! Record max of the surrounding ice thicknesses for correction in limupdate 
     125         ! Record max of the surrounding ice thicknesses for correction 
    118126         ! in case advection creates ice too thick. 
    119127         !--------------------------------------------------------------------- 
    … …  
    142150 
    143151         IF( zcfl > 0.5_wp .AND. lwp )   ncfl = ncfl + 1 
    144          IF( lwp ) THEN 
    145             IF( ncfl > 0 ) THEN    
    146                WRITE(cltmp,'(i6.1)') ncfl 
    147                CALL ctl_warn( 'lim_trp: ncfl= ', TRIM(cltmp), 'advective ice time-step using a split in sub-time-step ') 
    148             ELSE 
    149                WRITE(numout,*) 'lim_trp : CFL criterion for ice advection is always smaller than 1/2 ' 
    150             ENDIF 
    151          ENDIF 
     152!!         IF( lwp ) THEN 
     153!!            IF( ncfl > 0 ) THEN    
     154!!               WRITE(cltmp,'(i6.1)') ncfl 
     155!!               CALL ctl_warn( 'lim_trp: ncfl= ', TRIM(cltmp), 'advective ice time-step using a split in sub-time-step ') 
     156!!            ELSE 
     157!!            !  WRITE(numout,*) 'lim_trp : CFL criterion for ice advection is always smaller than 1/2 ' 
     158!!            ENDIF 
     159!!         ENDIF 
    152160 
    153161         !------------------------- 
    … …  
    228236                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   & 
    229237                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    230  
    231238                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0oi  (:,:,jl), sxage(:,:,jl),   &   !--- ice age      --- 
    232239                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    … …  
    345352!!gm & cr  
    346353 
     354         ! --- diags --- 
     355         DO jj = 1, jpj 
     356            DO ji = 1, jpi 
     357               diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) * r1_rdtice 
     358               diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) * r1_rdtice 
     359 
     360               diag_trp_vi (ji,jj) = SUM(   v_i(ji,jj,:) -  zviold(ji,jj,:) ) * r1_rdtice 
     361               diag_trp_vs (ji,jj) = SUM(   v_s(ji,jj,:) -  zvsold(ji,jj,:) ) * r1_rdtice 
     362               diag_trp_smv(ji,jj) = SUM( smv_i(ji,jj,:) - zsmvold(ji,jj,:) ) * r1_rdtice 
     363            END DO 
     364         END DO 
     365 
    347366         ! zap small areas 
    348367         CALL lim_var_zapsmall 
    349368 
    350369         !--- Thickness correction in case too high -------------------------------------------------------- 
    351          CALL lim_var_glo2eqv 
    352370         DO jl = 1, jpl 
    353371            DO jj = 1, jpj 
    … …  
    356374                  IF ( v_i(ji,jj,jl) > 0._wp ) THEN 
    357375 
     376                     rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) 
     377                     ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 
     378                     ht_s  (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 
     379                      
    358380                     zvi  = v_i  (ji,jj,jl) 
    359381                     zvs  = v_s  (ji,jj,jl) 
    … …  
    365387 
    366388                     IF ( ( zdv >  0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) .AND. zatold(ji,jj) < 0.80 ) .OR. & 
    367                         & ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN                                           
     389                        & ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN 
    368390 
    369391                        rswitch        = MAX( 0._wp, SIGN( 1._wp, zhimax(ji,jj,jl) - epsi20 ) ) 
    … …  
    405427         ENDIF 
    406428 
    407          ! --- diags --- 
    408          DO jj = 1, jpj 
    409             DO ji = 1, jpi 
    410                diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) * r1_rdtice 
    411                diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) * r1_rdtice 
    412  
    413                diag_trp_vi (ji,jj) = SUM(   v_i(ji,jj,:) -  zviold(ji,jj,:) ) * r1_rdtice 
    414                diag_trp_vs (ji,jj) = SUM(   v_s(ji,jj,:) -  zvsold(ji,jj,:) ) * r1_rdtice 
    415                diag_trp_smv(ji,jj) = SUM( smv_i(ji,jj,:) - zsmvold(ji,jj,:) ) * r1_rdtice 
    416             END DO 
    417          END DO 
    418  
    419429         ! --- agglomerate variables ----------------- 
    420430         vt_i (:,:) = 0._wp 
    … …  
    444454      ENDIF 
    445455 
    446       CALL lim_var_glo2eqv            ! equivalent variables, requested for rafting 
    447  
    448456      ! ------------------------------------------------- 
    449457      ! control prints 
    … …  
    451459      IF( ln_icectl )   CALL lim_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' ) 
    452460      ! 
    453       CALL wrk_dealloc( jpi,jpj,           zsm, zatold, zeiold, zesold ) 
    454       CALL wrk_dealloc( jpi,jpj,jpl,       z0ice, z0snw, z0ai, z0es , z0smi , z0oi ) 
    455       CALL wrk_dealloc( jpi,jpj,1,         z0opw ) 
    456       CALL wrk_dealloc( jpi,jpj,nlay_i+1,jpl, z0ei ) 
    457       CALL wrk_dealloc( jpi,jpj,jpl,       zviold, zvsold, zhimax, zsmvold ) 
     461      CALL wrk_dealloc( jpi,jpj,            zsm, zatold, zeiold, zesold ) 
     462      CALL wrk_dealloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi ) 
     463      CALL wrk_dealloc( jpi,jpj,1,          z0opw ) 
     464      CALL wrk_dealloc( jpi,jpj,nlay_i,jpl, z0ei ) 
     465      CALL wrk_dealloc( jpi,jpj,jpl,        zviold, zvsold, zhimax, zsmvold ) 
    458466      ! 
    459467      IF( nn_timing == 1 )  CALL timing_stop('limtrp') 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90

    • Property svn:keywords set to Id
    r5134 r5621  
    3939   !!---------------------------------------------------------------------- 
    4040   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    41    !! $Id: limupdate.F90 3294 2012-01-28 16:44:18Z rblod $ 
     41   !! $Id$ 
    4242   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4343   !!---------------------------------------------------------------------- 
    … …  
    6969      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    7070 
    71       CALL lim_var_glo2eqv 
    7271      !---------------------------------------------------- 
    7372      ! ice concentration should not exceed amax  
    … …  
    8281            DO ji = 1, jpi 
    8382               IF( at_i(ji,jj) > rn_amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN 
    84                   a_i(ji,jj,jl)  = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 
     83                  a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 
     84                  oa_i(ji,jj,jl) = oa_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 
    8585               ENDIF 
    8686            END DO 
    … …  
    8888      END DO 
    8989     
    90       !---------------------------------------------------- 
    91       ! Rebin categories with thickness out of bounds 
    92       !---------------------------------------------------- 
    93       IF ( jpl > 1 ) CALL lim_itd_th_reb(1, jpl) 
    94  
    95       !----------------- 
    96       ! zap small values 
    97       !----------------- 
    98       CALL lim_var_zapsmall 
    99  
    10090      !--------------------- 
    10191      ! Ice salinity bounds 
    … …  
    10696               DO ji = 1, jpi 
    10797                  zsal            = smv_i(ji,jj,jl) 
    108                   smv_i(ji,jj,jl) = sm_i(ji,jj,jl) * v_i(ji,jj,jl) 
    10998                  ! salinity stays in bounds 
    11099                  rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) ) 
    … …  
    117106      ENDIF 
    118107 
     108      !---------------------------------------------------- 
     109      ! Rebin categories with thickness out of bounds 
     110      !---------------------------------------------------- 
     111      IF ( jpl > 1 ) CALL lim_itd_th_reb(1, jpl) 
     112 
     113      !----------------- 
     114      ! zap small values 
     115      !----------------- 
     116      CALL lim_var_zapsmall 
     117 
     118      ! ------------------------------------------------- 
     119      ! Diagnostics 
     120      ! ------------------------------------------------- 
     121      DO jl  = 1, jpl 
     122         afx_dyn(:,:) = afx_dyn(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice 
     123      END DO 
     124 
     125      DO jj = 1, jpj 
     126         DO ji = 1, jpi             
     127            ! heat content variation (W.m-2) 
     128            diag_heat(ji,jj) = - ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  &  
     129               &                   SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    & 
     130               &                 ) * r1_rdtice 
     131            ! salt, volume 
     132            diag_smvi(ji,jj) = SUM( smv_i(ji,jj,:) - smv_i_b(ji,jj,:) ) * rhoic * r1_rdtice 
     133            diag_vice(ji,jj) = SUM( v_i  (ji,jj,:) - v_i_b  (ji,jj,:) ) * rhoic * r1_rdtice 
     134            diag_vsnw(ji,jj) = SUM( v_s  (ji,jj,:) - v_s_b  (ji,jj,:) ) * rhosn * r1_rdtice 
     135         END DO 
     136      END DO 
     137 
    119138      ! conservation test 
    120139      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    121  
    122       ! ------------------------------------------------- 
    123       ! Diagnostics 
    124       ! ------------------------------------------------- 
    125       DO jl  = 1, jpl 
    126          afx_dyn(:,:) = afx_dyn(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice 
    127       END DO 
    128  
    129       ! heat content variation (W.m-2) 
    130       DO jj = 1, jpj 
    131          DO ji = 1, jpi             
    132             diag_heat_dhc(ji,jj) = - ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  &  
    133                &                       SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    & 
    134                &                     ) * r1_rdtice    
    135          END DO 
    136       END DO 
    137140 
    138141      ! ------------------------------------------------- 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90

    • Property svn:keywords set to Id
    r5134 r5621  
    4141   !!---------------------------------------------------------------------- 
    4242   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    43    !! $Id: limupdate.F90 3294 2012-01-28 16:44:18Z rblod $ 
     43   !! $Id$ 
    4444   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4545   !!---------------------------------------------------------------------- 
    … …  
    7272      ! Constrain the thickness of the smallest category above himin 
    7373      !---------------------------------------------------------------------- 
    74       CALL lim_var_glo2eqv 
    7574      DO jj = 1, jpj  
    7675         DO ji = 1, jpi 
     76            rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,1) - epsi20 ) )   !0 if no ice and 1 if yes 
     77            ht_i(ji,jj,1) = v_i (ji,jj,1) / MAX( a_i(ji,jj,1) , epsi20 ) * rswitch 
    7778            IF( v_i(ji,jj,1) > 0._wp .AND. ht_i(ji,jj,1) < rn_himin ) THEN 
    78                a_i (ji,jj,1) = a_i(ji,jj,1) * ht_i(ji,jj,1) / rn_himin 
     79               a_i (ji,jj,1) = a_i (ji,jj,1) * ht_i(ji,jj,1) / rn_himin 
     80               oa_i(ji,jj,1) = oa_i(ji,jj,1) * ht_i(ji,jj,1) / rn_himin 
    7981            ENDIF 
    8082         END DO 
    … …  
    9395            DO ji = 1, jpi 
    9496               IF( at_i(ji,jj) > rn_amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN 
    95                   a_i(ji,jj,jl)  = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 
     97                  a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 
     98                  oa_i(ji,jj,jl) = oa_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 
    9699               ENDIF 
    97100            END DO 
    98101         END DO 
    99102      END DO 
    100  
    101       !---------------------------------------------------- 
    102       ! Rebin categories with thickness out of bounds 
    103       !---------------------------------------------------- 
    104       IF ( jpl > 1 ) CALL lim_itd_th_reb( 1, jpl ) 
    105  
    106       !----------------- 
    107       ! zap small values 
    108       !----------------- 
    109       CALL lim_var_zapsmall 
    110103 
    111104      !--------------------- 
    … …  
    117110               DO ji = 1, jpi 
    118111                  zsal            = smv_i(ji,jj,jl) 
    119                   smv_i(ji,jj,jl) = sm_i(ji,jj,jl) * v_i(ji,jj,jl) 
    120112                  ! salinity stays in bounds 
    121113                  rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) ) 
    … …  
    127119         END DO 
    128120      ENDIF 
     121 
     122      !---------------------------------------------------- 
     123      ! Rebin categories with thickness out of bounds 
     124      !---------------------------------------------------- 
     125      IF ( jpl > 1 ) CALL lim_itd_th_reb( 1, jpl ) 
     126 
     127      !----------------- 
     128      ! zap small values 
     129      !----------------- 
     130      CALL lim_var_zapsmall 
    129131 
    130132      !------------------------------------------------------------------------------ 
    … …  
    150152      v_ice(:,:) = v_ice(:,:) * vmask(:,:,1) 
    151153  
    152       ! for outputs 
    153       CALL lim_var_glo2eqv            ! equivalent variables (outputs) 
    154       CALL lim_var_agg(2)             ! aggregate ice thickness categories 
    155  
    156       ! conservation test 
    157       IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    158  
    159154      ! ------------------------------------------------- 
    160155      ! Diagnostics 
    161156      ! ------------------------------------------------- 
    162157      DO jl  = 1, jpl 
     158         oa_i(:,:,jl) = oa_i(:,:,jl) + a_i(:,:,jl) * rdt_ice / rday   ! ice natural aging 
    163159         afx_thd(:,:) = afx_thd(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice 
    164160      END DO 
    165161      afx_tot = afx_thd + afx_dyn 
    166162 
    167       ! heat content variation (W.m-2) 
    168163      DO jj = 1, jpj 
    169164         DO ji = 1, jpi             
    170             diag_heat_dhc(ji,jj) = diag_heat_dhc(ji,jj) -  & 
    171                &                   ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  &  
    172                &                     SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    & 
    173                &                   ) * r1_rdtice    
    174          END DO 
    175       END DO 
     165            ! heat content variation (W.m-2) 
     166            diag_heat(ji,jj) = diag_heat(ji,jj) -  & 
     167               &               ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  &  
     168               &                 SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    & 
     169               &               ) * r1_rdtice    
     170            ! salt, volume 
     171            diag_smvi(ji,jj) = diag_smvi(ji,jj) + SUM( smv_i(ji,jj,:) - smv_i_b(ji,jj,:) ) * rhoic * r1_rdtice 
     172            diag_vice(ji,jj) = diag_vice(ji,jj) + SUM( v_i  (ji,jj,:) - v_i_b  (ji,jj,:) ) * rhoic * r1_rdtice 
     173            diag_vsnw(ji,jj) = diag_vsnw(ji,jj) + SUM( v_s  (ji,jj,:) - v_s_b  (ji,jj,:) ) * rhosn * r1_rdtice 
     174         END DO 
     175      END DO 
     176 
     177      ! conservation test 
     178      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     179 
     180      ! necessary calls (at least for coupling) 
     181      CALL lim_var_glo2eqv 
     182      CALL lim_var_agg(2) 
    176183 
    177184      ! ------------------------------------------------- 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90

    r5134 r5621  
    124124               DO ji = 1, jpi 
    125125                  et_s(ji,jj)  = et_s(ji,jj)  + e_s(ji,jj,1,jl)                                           ! snow heat content 
    126                   rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )  
    127                   smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi10 ) * rswitch   ! ice salinity 
    128                   rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )  
    129                   ot_i(ji,jj)  = ot_i(ji,jj)  + oa_i(ji,jj,jl)  / MAX( at_i(ji,jj) , epsi10 ) * rswitch   ! ice age 
     126                  rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi20 ) )  
     127                  smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi20 ) * rswitch   ! ice salinity 
     128                  rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi20 ) )  
     129                  ot_i(ji,jj)  = ot_i(ji,jj)  + oa_i(ji,jj,jl)  / MAX( at_i(ji,jj) , epsi20 ) * rswitch   ! ice age 
    130130               END DO 
    131131            END DO 
    … …  
    161161         DO jj = 1, jpj 
    162162            DO ji = 1, jpi 
    163                rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )   !0 if no ice and 1 if yes 
    164                ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch 
    165                ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch 
    166                o_i(ji,jj,jl)  = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch 
     163               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )   !0 if no ice and 1 if yes 
     164               ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 
     165               ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 
     166               o_i(ji,jj,jl)  = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 
    167167            END DO 
    168168         END DO 
    … …  
    173173            DO jj = 1, jpj 
    174174               DO ji = 1, jpi 
    175                   rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) )   !0 if no ice and 1 if yes 
    176                   sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi10 ) * rswitch 
     175                  rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi20 ) )   !0 if no ice and 1 if yes 
     176                  sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * rswitch 
     177                  !                                      ! bounding salinity 
     178                  sm_i(ji,jj,jl) = MAX( sm_i(ji,jj,jl), rn_simin ) 
    177179               END DO 
    178180            END DO 
    … …  
    199201                  zdiscrim   =  SQRT( MAX(zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp) ) 
    200202                  t_i(ji,jj,jk,jl) = rt0 + rswitch *( - zbbb - zdiscrim ) / ( 2.0 *zaaa ) 
    201                   t_i(ji,jj,jk,jl) = MIN( rt0, MAX( 173.15_wp, t_i(ji,jj,jk,jl) ) )       ! 100-rt0 < t_i < rt0 
     203                  t_i(ji,jj,jk,jl) = MIN( ztmelts, MAX( rt0 - 100._wp, t_i(ji,jj,jk,jl) ) )  ! -100 < t_i < ztmelts 
    202204               END DO 
    203205            END DO 
    … …  
    219221                  ! 
    220222                  t_s(ji,jj,jk,jl) = rt0 + rswitch * ( - zfac1 * zq_s + zfac2 ) 
    221                   t_s(ji,jj,jk,jl) = MIN( rt0, MAX( 173.15, t_s(ji,jj,jk,jl) ) )     ! 100-rt0 < t_i < rt0 
     223                  t_s(ji,jj,jk,jl) = MIN( rt0, MAX( rt0 - 100._wp , t_s(ji,jj,jk,jl) ) )     ! -100 < t_s < rt0 
    222224               END DO 
    223225            END DO 
    … …  
    228230      ! Mean temperature 
    229231      !------------------- 
     232      vt_i (:,:) = 0._wp 
     233      DO jl = 1, jpl 
     234         vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 
     235      END DO 
     236 
    230237      tm_i(:,:) = 0._wp 
    231238      DO jl = 1, jpl 
    … …  
    234241               DO ji = 1, jpi 
    235242                  rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 
    236                   tm_i(ji,jj) = tm_i(ji,jj) + rswitch * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl)   & 
    237                      &                      / (  REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 )  ) 
    238                END DO 
    239             END DO 
    240          END DO 
    241       END DO 
     243                  tm_i(ji,jj) = tm_i(ji,jj) + r1_nlay_i * rswitch * ( t_i(ji,jj,jk,jl) - rt0 ) * v_i(ji,jj,jl)  & 
     244                     &            / MAX( vt_i(ji,jj) , epsi10 ) 
     245               END DO 
     246            END DO 
     247         END DO 
     248      END DO 
     249      tm_i = tm_i + rt0 
    242250      ! 
    243251   END SUBROUTINE lim_var_glo2eqv 
    … …  
    258266      v_s(:,:,:)   = ht_s(:,:,:) * a_i(:,:,:) 
    259267      smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:) 
    260       oa_i (:,:,:) = o_i (:,:,:) * a_i(:,:,:) 
    261268      ! 
    262269   END SUBROUTINE lim_var_eqv2glo 
    … …  
    305312            DO jj = 1, jpj 
    306313               DO ji = 1, jpi 
    307                   z_slope_s(ji,jj,jl) = 2._wp * sm_i(ji,jj,jl) / MAX( epsi10 , ht_i(ji,jj,jl) ) 
     314                  rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i(ji,jj,jl) - epsi20 ) ) 
     315                  z_slope_s(ji,jj,jl) = rswitch * 2._wp * sm_i(ji,jj,jl) / MAX( epsi20 , ht_i(ji,jj,jl) ) 
    308316               END DO 
    309317            END DO 
    … …  
    339347                     !                                      ! weighting the profile 
    340348                     s_i(ji,jj,jk,jl) = zalpha(ji,jj,jl) * zs_zero + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl) 
     349                     !                                      ! bounding salinity 
     350                     s_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( s_i(ji,jj,jk,jl), rn_simin ) ) 
    341351                  END DO 
    342352               END DO 
    … …  
    379389 
    380390      ! Mean sea ice temperature 
     391      vt_i (:,:) = 0._wp 
     392      DO jl = 1, jpl 
     393         vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 
     394      END DO 
     395 
    381396      tm_i(:,:) = 0._wp 
    382397      DO jl = 1, jpl 
    … …  
    385400               DO ji = 1, jpi 
    386401                  rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 
    387                   tm_i(ji,jj) = tm_i(ji,jj) + rswitch * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl)   & 
    388                      &                      * r1_nlay_i / MAX( vt_i(ji,jj) , epsi10 ) 
    389                END DO 
    390             END DO 
    391          END DO 
    392       END DO 
     402                  tm_i(ji,jj) = tm_i(ji,jj) + r1_nlay_i * rswitch * ( t_i(ji,jj,jk,jl) - rt0 ) * v_i(ji,jj,jl)  & 
     403                     &            / MAX( vt_i(ji,jj) , epsi10 ) 
     404               END DO 
     405            END DO 
     406         END DO 
     407      END DO 
     408      tm_i = tm_i + rt0 
    393409 
    394410   END SUBROUTINE lim_var_icetm 
    … …  
    409425      !!------------------------------------------------------------------ 
    410426      ! 
     427      vt_i (:,:) = 0._wp 
     428      DO jl = 1, jpl 
     429         vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 
     430      END DO 
     431 
    411432      bv_i(:,:) = 0._wp 
    412433      DO jl = 1, jpl 
    … …  
    417438                  zbvi  = - rswitch * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 )   & 
    418439                     &                   * v_i(ji,jj,jl) * r1_nlay_i 
    419                   rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  ) 
    420                   bv_i(ji,jj) = bv_i(ji,jj) + rswitch * zbvi  / MAX( vt_i(ji,jj) , epsi10 ) 
     440                  rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi20 ) )  ) 
     441                  bv_i(ji,jj) = bv_i(ji,jj) + rswitch * zbvi  / MAX( vt_i(ji,jj) , epsi20 ) 
    421442               END DO 
    422443            END DO 
    … …  
    460481         ! 
    461482         DO ji = kideb, kiut          ! Slope of the linear profile zs_zero 
    462             z_slope_s(ji) = 2._wp * sm_i_1d(ji) / MAX( epsi10 , ht_i_1d(ji) ) 
     483            rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) ) 
     484            z_slope_s(ji) = rswitch * 2._wp * sm_i_1d(ji) / MAX( epsi20 , ht_i_1d(ji) ) 
    463485         END DO 
    464486 
    … …  
    484506               ! weighting the profile 
    485507               s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji) 
     508               ! bounding salinity 
     509               s_i_1d(ji,jk) = MIN( rn_simax, MAX( s_i_1d(ji,jk), rn_simin ) ) 
    486510            END DO  
    487511         END DO  
    … …  
    537561                  rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) 
    538562                  rswitch          = MAX( 0._wp , SIGN( 1._wp, at_i(ji,jj  ) - epsi10 ) ) * rswitch 
     563                  rswitch          = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch 
     564                  rswitch          = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch  & 
     565                     &                                       / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch 
    539566                  zei              = e_i(ji,jj,jk,jl) 
    540567                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch 
    … …  
    550577               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) 
    551578               rswitch = MAX( 0._wp , SIGN( 1._wp, at_i(ji,jj  ) - epsi10 ) ) * rswitch 
    552                 
     579               rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch 
     580               rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch  & 
     581                  &                              / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch 
    553582               zsal = smv_i(ji,jj,  jl) 
    554583               zvi  = v_i  (ji,jj,  jl) 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90

    r5123 r5621  
    6060      REAL(wp) ::  z1_365 
    6161      REAL(wp) ::  ztmp 
    62       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zoi, zei 
     62      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zoi, zei, zt_i, zt_s 
    6363      REAL(wp), POINTER, DIMENSION(:,:)   ::  z2d, z2da, z2db, zswi    ! 2D workspace 
    6464      !!------------------------------------------------------------------- 
    … …  
    6666      IF( nn_timing == 1 )  CALL timing_start('limwri') 
    6767 
    68       CALL wrk_alloc( jpi, jpj, jpl, zoi, zei ) 
     68      CALL wrk_alloc( jpi, jpj, jpl, zoi, zei, zt_i, zt_s ) 
    6969      CALL wrk_alloc( jpi, jpj     , z2d, z2da, z2db, zswi ) 
    7070 
    … …  
    7272      ! Mean category values 
    7373      !----------------------------- 
     74      z1_365 = 1._wp / 365._wp 
    7475 
    7576      CALL lim_var_icetm      ! mean sea ice temperature 
    … …  
    112113         CALL lbc_lnk( z2da, 'T', -1. ) 
    113114         CALL lbc_lnk( z2db, 'T', -1. ) 
    114          CALL iom_put( "uice_ipa"     , z2da                )       ! ice velocity u component 
    115          CALL iom_put( "vice_ipa"     , z2db                )       ! ice velocity v component 
     115         CALL iom_put( "uice_ipa"     , z2da             )       ! ice velocity u component 
     116         CALL iom_put( "vice_ipa"     , z2db             )       ! ice velocity v component 
    116117         DO jj = 1, jpj                                  
    117118            DO ji = 1, jpi 
    … …  
    119120            END DO 
    120121         END DO 
    121          CALL iom_put( "icevel"       , z2d                 )       ! ice velocity module 
     122         CALL iom_put( "icevel"       , z2d              )       ! ice velocity module 
    122123      ENDIF 
    123124      ! 
    … …  
    127128            DO jj = 1, jpj 
    128129               DO ji = 1, jpi 
    129                   z2d(ji,jj) = z2d(ji,jj) + zswi(ji,jj) * oa_i(ji,jj,jl) 
     130                  rswitch    = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.1 ) ) 
     131                  z2d(ji,jj) = z2d(ji,jj) + rswitch * oa_i(ji,jj,jl) / MAX( at_i(ji,jj), 0.1 ) 
    130132               END DO 
    131133            END DO 
    132134         END DO 
    133          z1_365 = 1._wp / 365._wp 
    134          CALL iom_put( "miceage"     , z2d * z1_365         )        ! mean ice age 
     135         CALL iom_put( "miceage"     , z2d * z1_365      )        ! mean ice age 
    135136      ENDIF 
    136137 
    … …  
    141142            END DO 
    142143         END DO 
    143          CALL iom_put( "micet"       , z2d                  )        ! mean ice temperature 
     144         CALL iom_put( "micet"       , z2d               )        ! mean ice temperature 
    144145      ENDIF 
    145146      ! 
    … …  
    153154            END DO 
    154155         END DO 
    155          CALL iom_put( "icest"       , z2d                 )        ! ice surface temperature 
     156         CALL iom_put( "icest"       , z2d              )        ! ice surface temperature 
    156157      ENDIF 
    157158 
    … …  
    163164            END DO 
    164165         END DO 
    165          CALL iom_put( "icecolf"     , z2d                 )        ! frazil ice collection thickness 
     166         CALL iom_put( "icecolf"     , z2d              )        ! frazil ice collection thickness 
    166167      ENDIF 
    167168 
    … …  
    175176      CALL iom_put( "utau_ice"    , utau_ice            )        ! wind stress over ice along i-axis at I-point 
    176177      CALL iom_put( "vtau_ice"    , vtau_ice            )        ! wind stress over ice along j-axis at I-point 
    177       CALL iom_put( "snowpre"     , sprecip             )        ! snow precipitation  
     178      CALL iom_put( "snowpre"     , sprecip * 86400.    )        ! snow precipitation  
    178179      CALL iom_put( "micesalt"    , smt_i               )        ! mean ice salinity 
    179180 
    … …  
    231232      CALL iom_put ('hfxdif'     , hfx_dif(:,:)         )   !   
    232233      CALL iom_put ('hfxopw'     , hfx_opw(:,:)         )   !   
    233       CALL iom_put ('hfxtur'     , fhtur(:,:) * at_i(:,:) ) ! turbulent heat flux at ice base  
    234       CALL iom_put ('hfxdhc'     , diag_heat_dhc(:,:)   )   ! Heat content variation in snow and ice  
     234      CALL iom_put ('hfxtur'     , fhtur(:,:) * SUM(a_i_b(:,:,:), dim=3) ) ! turbulent heat flux at ice base  
     235      CALL iom_put ('hfxdhc'     , diag_heat(:,:)       )   ! Heat content variation in snow and ice  
    235236      CALL iom_put ('hfxspr'     , hfx_spr(:,:)         )   ! Heat content of snow precip  
    236237       
    … …  
    243244      CALL iom_put( "salinity_cat"     , sm_i        )        ! salinity for categories 
    244245 
     246      ! ice temperature 
     247      IF ( iom_use( "icetemp_cat" ) ) THEN  
     248         zt_i(:,:,:) = SUM( t_i(:,:,:,:), dim=3 ) * r1_nlay_i 
     249         CALL iom_put( "icetemp_cat"   , zt_i - rt0  ) 
     250      ENDIF 
     251       
     252      ! snow temperature 
     253      IF ( iom_use( "snwtemp_cat" ) ) THEN  
     254         zt_s(:,:,:) = SUM( t_s(:,:,:,:), dim=3 ) * r1_nlay_s 
     255         CALL iom_put( "snwtemp_cat"   , zt_s - rt0  ) 
     256      ENDIF 
     257 
    245258      ! Compute ice age 
    246259      IF ( iom_use( "iceage_cat" ) ) THEN  
    … …  
    248261            DO jj = 1, jpj 
    249262               DO ji = 1, jpi 
    250                   rswitch = MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - epsi06 ) ) 
    251                   zoi(ji,jj,jl) = oa_i(ji,jj,jl)  / MAX( a_i(ji,jj,jl) , epsi06 ) * rswitch 
     263                  rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.1 ) ) 
     264                  rswitch = rswitch * MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - 0.1 ) ) 
     265                  zoi(ji,jj,jl) = oa_i(ji,jj,jl)  / MAX( a_i(ji,jj,jl) , 0.1 ) * rswitch 
    252266               END DO 
    253267            END DO 
    254268         END DO 
    255          CALL iom_put( "iceage_cat"     , zoi        )        ! ice age for categories 
     269         CALL iom_put( "iceage_cat"   , zoi * z1_365 )        ! ice age for categories 
    256270      ENDIF 
    257271 
    … …  
    264278                  DO ji = 1, jpi 
    265279                     rswitch = MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - epsi06 ) ) 
    266                      zei(ji,jj,jl) = zei(ji,jj,jl) + 100.0* & 
     280                     zei(ji,jj,jl) = zei(ji,jj,jl) + 100.0 * & 
    267281                        ( - tmut * s_i(ji,jj,jk,jl) / MIN( ( t_i(ji,jj,jk,jl) - rt0 ), - epsi06 ) ) * & 
    268282                        rswitch * r1_nlay_i 
    … …  
    271285            END DO 
    272286         END DO 
    273          CALL iom_put( "brinevol_cat"     , zei         )        ! brine volume for categories 
     287         CALL iom_put( "brinevol_cat"     , zei      )        ! brine volume for categories 
    274288      ENDIF 
    275289 
    … …  
    278292      !     not yet implemented 
    279293       
    280       CALL wrk_dealloc( jpi, jpj, jpl, zoi, zei ) 
     294      CALL wrk_dealloc( jpi, jpj, jpl, zoi, zei, zt_i, zt_s ) 
    281295      CALL wrk_dealloc( jpi, jpj     , z2d, zswi, z2da, z2db ) 
    282296 

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90

    r5146 r5621  
    2020   !                               !!! ** ice-thermo namelist (namicethd) ** 
    2121   REAL(wp), PUBLIC ::   rn_himin    !: minimum ice thickness 
    22    REAL(wp), PUBLIC ::   parsub      !: switch for snow sublimation or not 
    2322   REAL(wp), PUBLIC ::   rn_maxfrazb !: maximum portion of frazil ice collecting at the ice bottom 
    2423   REAL(wp), PUBLIC ::   rn_vfrazb   !: threshold drift speed for collection of bottom frazil ice 
    … …  
    9089   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fhld_1d       !: <==> the 2D  fhld 
    9190   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dqns_ice_1d   !: <==> the 2D  dqns_ice 
    92    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qla_ice_1d    !: <==> the 2D  qla_ice 
    93    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dqla_ice_1d   !: <==> the 2D  dqla_ice 
    94    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   tatm_ice_1d   !: <==> the 2D  tatm_ice 
     91   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   evap_ice_1d   !: <==> the 2D  evap_ice 
     92   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qprec_ice_1d  !: <==> the 2D  qprec_ice 
    9593   !                                                     ! to reintegrate longwave flux inside the ice thermodynamics 
    9694   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   i0            !: fraction of radiation transmitted to the ice 
    … …  
    140138      !!---------------------------------------------------------------------! 
    141139 
    142       ALLOCATE( npb      (jpij) , nplm        (jpij) , npac     (jpij),   & 
    143          !                                                                  ! 
    144          &      qlead_1d (jpij) , ftr_ice_1d  (jpij) ,     & 
    145          &      qsr_ice_1d (jpij) ,     & 
    146          &      fr1_i0_1d(jpij) , fr2_i0_1d(jpij) , qns_ice_1d(jpij) ,     & 
     140      ALLOCATE( npb      (jpij) , nplm      (jpij) , npac       (jpij) ,   & 
     141         &      qlead_1d (jpij) , ftr_ice_1d(jpij) , qsr_ice_1d (jpij) ,   & 
     142         &      fr1_i0_1d(jpij) , fr2_i0_1d (jpij) , qns_ice_1d(jpij)  ,   & 
    147143         &      t_bo_1d   (jpij) ,                                         & 
    148144         &      hfx_sum_1d(jpij) , hfx_bom_1d(jpij) ,hfx_bog_1d(jpij) ,    &  
    … …  
    152148         &      hfx_res_1d(jpij) , hfx_err_rem_1d(jpij) , hfx_err_dif_1d(jpij) , STAT=ierr(1) ) 
    153149      ! 
    154       ALLOCATE( sprecip_1d (jpij) , frld_1d    (jpij) , at_i_1d     (jpij) ,     & 
    155          &      fhtur_1d   (jpij) , wfx_snw_1d (jpij) , wfx_spr_1d (jpij) ,     & 
    156          &      fhld_1d    (jpij) , wfx_sub_1d (jpij) , wfx_bog_1d(jpij) , wfx_bom_1d(jpij) , & 
    157          &      wfx_sum_1d(jpij)  , wfx_sni_1d (jpij) , wfx_opw_1d (jpij) ,  wfx_res_1d (jpij) ,  & 
    158          &      dqns_ice_1d(jpij) , qla_ice_1d (jpij) , dqla_ice_1d(jpij) ,     & 
    159          &      tatm_ice_1d(jpij) ,      &    
    160          &      i0         (jpij) ,     &   
    161          &      sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) ,sfx_sum_1d (jpij) ,   & 
    162          &      sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , & 
    163          &      dsm_i_fl_1d(jpij) , dsm_i_gd_1d(jpij) , dsm_i_se_1d(jpij) ,     &      
     150      ALLOCATE( sprecip_1d (jpij) , frld_1d    (jpij) , at_i_1d    (jpij) ,                     & 
     151         &      fhtur_1d   (jpij) , wfx_snw_1d (jpij) , wfx_spr_1d (jpij) ,                     & 
     152         &      fhld_1d    (jpij) , wfx_sub_1d (jpij) , wfx_bog_1d (jpij) , wfx_bom_1d(jpij) ,  & 
     153         &      wfx_sum_1d(jpij)  , wfx_sni_1d (jpij) , wfx_opw_1d (jpij) , wfx_res_1d(jpij) ,  & 
     154         &      dqns_ice_1d(jpij) , evap_ice_1d (jpij),                                         & 
     155         &      qprec_ice_1d(jpij), i0         (jpij) ,                                         &   
     156         &      sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) , sfx_sum_1d (jpij),  & 
     157         &      sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) ,                     & 
     158         &      dsm_i_fl_1d(jpij) , dsm_i_gd_1d(jpij) , dsm_i_se_1d(jpij) ,                     &      
    164159         &      dsm_i_si_1d(jpij) , hicol_1d    (jpij)                     , STAT=ierr(2) ) 
    165160      ! 
    166       ALLOCATE( t_su_1d    (jpij) , a_i_1d    (jpij) , ht_i_1d   (jpij) ,    &    
    167          &      ht_s_1d    (jpij) , fc_su    (jpij) , fc_bo_i  (jpij) ,    &     
     161      ALLOCATE( t_su_1d   (jpij) , a_i_1d   (jpij) , ht_i_1d  (jpij) ,    &    
     162         &      ht_s_1d   (jpij) , fc_su    (jpij) , fc_bo_i  (jpij) ,    &     
    168163         &      dh_s_tot  (jpij) , dh_i_surf(jpij) , dh_i_bott(jpij) ,    &     
    169          &      dh_snowice(jpij) ,  & 
    170          &      sm_i_1d   (jpij) , s_i_new  (jpij) ,    & 
    171          &      t_s_1d(jpij,nlay_s),                                       & 
    172          &      t_i_1d(jpij,nlay_i+1), s_i_1d(jpij,nlay_i+1)                ,     &             
    173          &      q_i_1d(jpij,nlay_i+1), q_s_1d(jpij,nlay_i+1)                ,     & 
     164         &      dh_snowice(jpij) , sm_i_1d  (jpij) , s_i_new  (jpij) ,    & 
     165         &      t_s_1d(jpij,nlay_s) , t_i_1d(jpij,nlay_i) , s_i_1d(jpij,nlay_i) ,  &             
     166         &      q_i_1d(jpij,nlay_i+1) , q_s_1d(jpij,nlay_s) ,                        & 
    174167         &      qh_i_old(jpij,0:nlay_i+1), h_i_old(jpij,0:nlay_i+1) , STAT=ierr(3)) 
    175168      ! 
Note: See TracChangeset for help on using the changeset viewer.