Changeset 5064


Ignore:
Timestamp:
2015-02-05T18:54:24+01:00 (6 years ago)
Author:
clem
Message:

LIM3 now includes specific physics to run with only 1 sea ice category (i.e. LIM2 type)

Location:
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO
Files:
22 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/dom_ice.F90

    r5059 r5064  
    2121 
    2222   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fcor       !: coriolis coefficient 
    23    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   area       !: surface of grid cell  
    2423   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tms, tmi   !: temperature mask, mask for stress 
    2524   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmu, tmv   !: mask at u and v velocity points 
     
    4241      !!------------------------------------------------------------------- 
    4342      ! 
    44       ALLOCATE( fcor(jpi,jpj)   , area(jpi,jpj) ,      & 
     43      ALLOCATE( fcor(jpi,jpj)   ,                      & 
    4544         &      tms   (jpi,jpj) , tmi (jpi,jpj) ,      & 
    4645         &      tmu   (jpi,jpj) , tmv (jpi,jpj) ,      & 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r5059 r5064  
    109109   !! smv_i       |      -      |    Sea ice salt content         | ppt.m | 
    110110   !! oa_i        !      -      !    Sea ice areal age content    | day   | 
    111    !! e_i         !      -      !    Ice enthalpy                 | 10^9 J|  
     111   !! e_i         !      -      !    Ice enthalpy                 | J/m2  |  
    112112   !!      -      ! q_i_1d      !    Ice enthalpy per unit vol.   | J/m3  |  
    113    !! e_s         !      -      !    Snow enthalpy                | 10^9 J|  
     113   !! e_s         !      -      !    Snow enthalpy                | J/m2  |  
    114114   !!      -      ! q_s_1d      !    Snow enthalpy per unit vol.  | J/m3  |  
    115115   !!                                                                     | 
     
    224224 
    225225   !                                     !!** define some parameters  
    226    REAL(wp), PUBLIC, PARAMETER ::   unit_fac = 1.e+09_wp  !: conversion factor for ice / snow enthalpy 
    227226   REAL(wp), PUBLIC, PARAMETER ::   epsi06   = 1.e-06_wp  !: small number  
    228227   REAL(wp), PUBLIC, PARAMETER ::   epsi10   = 1.e-10_wp  !: small number  
     
    331330       
    332331   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   t_i        !: ice temperatures          [K] 
    333    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i        !: ice thermal contents [Giga J] 
     332   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i        !: ice thermal contents    [J/m2] 
    334333   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   s_i        !: ice salinities          [PSU] 
    335334 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90

    r4990 r5064  
    9797 
    9898      !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                      
    99       psm (:,:)  = MAX( pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 ) 
     99      psm (:,:)  = MAX( pcrh * e12t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 ) 
    100100 
    101101      !  Calculate fluxes and moments between boxes i<-->i+1               
     
    282282 
    283283      !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 
    284       psm(:,:)  = MAX(  pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20  ) 
     284      psm(:,:)  = MAX(  pcrh * e12t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20  ) 
    285285 
    286286      !  Calculate fluxes and moments between boxes j<-->j+1               
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90

    r5059 r5064  
    167167      REAL(wp)                        :: zvi,   zsmv,   zei,   zfs,   zfw,   zft 
    168168      REAL(wp)                        :: zvmin, zamin, zamax  
     169      REAL(wp)                        :: zconv 
     170 
     171      zconv = 1.e-9 
    169172 
    170173      IF( icount == 0 ) THEN 
    171174 
    172          zvi_b  = glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) 
    173          zsmv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
    174          zei_b  = glob_sum( SUM(   e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) 
     175         zvi_b  = glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * e12t(:,:) * tms(:,:) ) 
     176         zsmv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tms(:,:) ) 
     177         zei_b  = glob_sum( SUM( SUM(   e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) * e12t(:,:) * tms * zconv +  & 
     178            &               SUM( SUM(   e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) * e12t(:,:) * tms * zconv ) 
    175179         zfw_b  = glob_sum( - ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  & 
    176180            &                   wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:)    & 
    177             &             ) * area(:,:) * tms(:,:) ) 
     181            &             ) * e12t(:,:) * tms(:,:) ) 
    178182         zfs_b  = glob_sum(   ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
    179183            &                   sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:)                                  & 
    180             &                 ) * area(:,:) * tms(:,:) ) 
     184            &                 ) * e12t(:,:) * tms(:,:) ) 
    181185         zft_b  = glob_sum(   ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  &  
    182186            &                 - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   & 
    183             &                  ) * area(:,:) / unit_fac * tms(:,:) ) 
     187            &                  ) * e12t(:,:) * tms(:,:) * zconv ) 
    184188 
    185189      ELSEIF( icount == 1 ) THEN 
     
    187191         zfs  = glob_sum(   ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
    188192            &                 sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:)                                  &  
    189             &                ) * area(:,:) * tms(:,:) ) - zfs_b 
     193            &                ) * e12t(:,:) * tms(:,:) ) - zfs_b 
    190194         zfw  = glob_sum( - ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  & 
    191195            &                 wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:)    & 
    192             &                ) * area(:,:) * tms(:,:) ) - zfw_b 
     196            &                ) * e12t(:,:) * tms(:,:) ) - zfw_b 
    193197         zft  = glob_sum(   ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  &  
    194198            &               - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   & 
    195             &                ) * area(:,:) / unit_fac * tms(:,:) ) - zft_b 
     199            &                ) * e12t(:,:) * tms(:,:) * zconv ) - zft_b 
    196200  
    197          zvi  = ( glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zvi_b ) * r1_rdtice - zfw  
    198          zsmv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zsmv_b ) * r1_rdtice + ( zfs / rhoic ) 
    199          zei  =   glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zei_b * r1_rdtice + zft 
     201         zvi  = ( glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * e12t(:,:) * tms(:,:) ) - zvi_b ) * r1_rdtice - zfw  
     202         zsmv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tms(:,:) ) - zsmv_b ) * r1_rdtice + ( zfs / rhoic ) 
     203         zei  =   glob_sum( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) * e12t(:,:) * tms * zconv +   & 
     204            &               SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) * e12t(:,:) * tms * zconv ) * r1_rdtice - zei_b * r1_rdtice + zft 
    200205 
    201206         zvmin = glob_min(v_i) 
     
    206211            IF ( ABS( zvi    ) >  1.e-4 ) WRITE(numout,*) 'violation volume [kg/day]     (',cd_routine,') = ',(zvi * rday) 
    207212            IF ( ABS( zsmv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (',cd_routine,') = ',(zsmv * rday) 
    208             IF ( ABS( zei    ) >  1.    ) WRITE(numout,*) 'violation enthalpy [1e9 J]    (',cd_routine,') = ',(zei) 
    209             IF ( zvmin <  1.e-10        ) WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',(zvmin) 
    210             IF( cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' .AND. zamax > amax+1.e-10 ) THEN 
     213            IF ( ABS( zei    ) >  1.e-2 ) WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',(zei) 
     214            IF ( zvmin <  -epsi10       ) WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',(zvmin) 
     215            IF( cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' .AND. zamax > amax+epsi10 ) THEN 
    211216                                          WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax 
    212217            ENDIF 
    213             IF ( zamin <  1.e-10        ) WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin 
     218            IF ( zamin <  -epsi10       ) WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin 
    214219         ENDIF 
    215220 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90

    r5059 r5064  
    306306               WRITE(numout,*) ' - Cell values ' 
    307307               WRITE(numout,*) '   ~~~~~~~~~~~ ' 
    308                WRITE(numout,*) ' cell area     : ', area(ji,jj) 
     308               WRITE(numout,*) ' cell area     : ', e12t(ji,jj) 
    309309               WRITE(numout,*) ' at_i          : ', at_i(ji,jj)        
    310310               WRITE(numout,*) ' vt_i          : ', vt_i(ji,jj)        
     
    317317                  WRITE(numout,*) ' v_i           : ', v_i(ji,jj,jl) 
    318318                  WRITE(numout,*) ' v_s           : ', v_s(ji,jj,jl) 
    319                   WRITE(numout,*) ' e_s           : ', e_s(ji,jj,1,jl)/1.0e9 
    320                   WRITE(numout,*) ' e_i           : ', e_i(ji,jj,1:nlay_i,jl)/1.0e9 
     319                  WRITE(numout,*) ' e_s           : ', e_s(ji,jj,1,jl) 
     320                  WRITE(numout,*) ' e_i           : ', e_i(ji,jj,1:nlay_i,jl) 
    321321                  WRITE(numout,*) ' t_su          : ', t_su(ji,jj,jl) 
    322322                  WRITE(numout,*) ' t_snow        : ', t_s(ji,jj,1,jl) 
     
    350350               WRITE(numout,*) ' - Cell values ' 
    351351               WRITE(numout,*) '   ~~~~~~~~~~~ ' 
    352                WRITE(numout,*) ' cell area     : ', area(ji,jj) 
     352               WRITE(numout,*) ' cell area     : ', e12t(ji,jj) 
    353353               WRITE(numout,*) ' at_i          : ', at_i(ji,jj)        
    354354               WRITE(numout,*) ' vt_i          : ', vt_i(ji,jj)        
     
    372372                  WRITE(numout,*) ' v_i        : ', v_i(ji,jj,jl)              , ' v_i_b      : ', v_i_b(ji,jj,jl)    
    373373                  WRITE(numout,*) ' v_s        : ', v_s(ji,jj,jl)              , ' v_s_b      : ', v_s_b(ji,jj,jl)   
    374                   WRITE(numout,*) ' e_i1       : ', e_i(ji,jj,1,jl)/1.0e9      , ' ei1        : ', e_i_b(ji,jj,1,jl)/1.0e9  
    375                   WRITE(numout,*) ' e_i2       : ', e_i(ji,jj,2,jl)/1.0e9      , ' ei2_b      : ', e_i_b(ji,jj,2,jl)/1.0e9   
     374                  WRITE(numout,*) ' e_i1       : ', e_i(ji,jj,1,jl)            , ' ei1        : ', e_i_b(ji,jj,1,jl)  
     375                  WRITE(numout,*) ' e_i2       : ', e_i(ji,jj,2,jl)            , ' ei2_b      : ', e_i_b(ji,jj,2,jl)   
    376376                  WRITE(numout,*) ' e_snow     : ', e_s(ji,jj,1,jl)            , ' e_snow_b   : ', e_s_b(ji,jj,1,jl)  
    377377                  WRITE(numout,*) ' smv_i      : ', smv_i(ji,jj,jl)            , ' smv_i_b    : ', smv_i_b(ji,jj,jl)    
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90

    r5055 r5064  
    7171 
    7272      ! 1/area 
    73       z1_area = 1._wp / MAX( glob_sum( area(:,:) * tms(:,:) ), epsi06 ) 
    74  
    75       rswitch = MAX( 0._wp , SIGN( 1._wp , glob_sum( area(:,:) * tms(:,:) ) - epsi06 ) ) 
     73      z1_area = 1._wp / MAX( glob_sum( e12t(:,:) * tms(:,:) ), epsi06 ) 
     74 
     75      rswitch = MAX( 0._wp , SIGN( 1._wp , glob_sum( e12t(:,:) * tms(:,:) ) - epsi06 ) ) 
    7676      ! ----------------------- ! 
    7777      ! 1 -  Content variations ! 
    7878      ! ----------------------- ! 
    79       zbg_ivo = glob_sum( vt_i(:,:) * area(:,:) * tms(:,:) ) ! volume ice  
    80       zbg_svo = glob_sum( vt_s(:,:) * area(:,:) * tms(:,:) ) ! volume snow 
    81       zbg_are = glob_sum( at_i(:,:) * area(:,:) * tms(:,:) ) ! area 
    82       zbg_sal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )       ! mean salt content 
    83       zbg_tem = glob_sum( ( tm_i(:,:) - rtt ) * vt_i(:,:) * area(:,:) * tms(:,:) )  ! mean temp content 
    84  
    85       !zbg_ihc = glob_sum( et_i(:,:) * area(:,:) * tms(:,:) ) / MAX( zbg_ivo,epsi06 ) ! ice heat content 
    86       !zbg_shc = glob_sum( et_s(:,:) * area(:,:) * tms(:,:) ) / MAX( zbg_svo,epsi06 ) ! snow heat content 
     79      zbg_ivo = glob_sum( vt_i(:,:) * e12t(:,:) * tms(:,:) ) ! volume ice  
     80      zbg_svo = glob_sum( vt_s(:,:) * e12t(:,:) * tms(:,:) ) ! volume snow 
     81      zbg_are = glob_sum( at_i(:,:) * e12t(:,:) * tms(:,:) ) ! area 
     82      zbg_sal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tms(:,:) )       ! mean salt content 
     83      zbg_tem = glob_sum( ( tm_i(:,:) - rtt ) * vt_i(:,:) * e12t(:,:) * tms(:,:) )  ! mean temp content 
     84 
     85      !zbg_ihc = glob_sum( et_i(:,:) * e12t(:,:) * tms(:,:) ) / MAX( zbg_ivo,epsi06 ) ! ice heat content 
     86      !zbg_shc = glob_sum( et_s(:,:) * e12t(:,:) * tms(:,:) ) / MAX( zbg_svo,epsi06 ) ! snow heat content 
    8787 
    8888      ! Volume 
    8989      ztmp = rswitch * z1_area * r1_rau0 * rday 
    90       zbg_vfx     = ztmp * glob_sum(     emp(:,:) * area(:,:) * tms(:,:) ) 
    91       zbg_vfx_bog = ztmp * glob_sum( wfx_bog(:,:) * area(:,:) * tms(:,:) ) 
    92       zbg_vfx_opw = ztmp * glob_sum( wfx_opw(:,:) * area(:,:) * tms(:,:) ) 
    93       zbg_vfx_sni = ztmp * glob_sum( wfx_sni(:,:) * area(:,:) * tms(:,:) ) 
    94       zbg_vfx_dyn = ztmp * glob_sum( wfx_dyn(:,:) * area(:,:) * tms(:,:) ) 
    95       zbg_vfx_bom = ztmp * glob_sum( wfx_bom(:,:) * area(:,:) * tms(:,:) ) 
    96       zbg_vfx_sum = ztmp * glob_sum( wfx_sum(:,:) * area(:,:) * tms(:,:) ) 
    97       zbg_vfx_res = ztmp * glob_sum( wfx_res(:,:) * area(:,:) * tms(:,:) ) 
    98       zbg_vfx_spr = ztmp * glob_sum( wfx_spr(:,:) * area(:,:) * tms(:,:) ) 
    99       zbg_vfx_snw = ztmp * glob_sum( wfx_snw(:,:) * area(:,:) * tms(:,:) ) 
    100       zbg_vfx_sub = ztmp * glob_sum( wfx_sub(:,:) * area(:,:) * tms(:,:) ) 
     90      zbg_vfx     = ztmp * glob_sum(     emp(:,:) * e12t(:,:) * tms(:,:) ) 
     91      zbg_vfx_bog = ztmp * glob_sum( wfx_bog(:,:) * e12t(:,:) * tms(:,:) ) 
     92      zbg_vfx_opw = ztmp * glob_sum( wfx_opw(:,:) * e12t(:,:) * tms(:,:) ) 
     93      zbg_vfx_sni = ztmp * glob_sum( wfx_sni(:,:) * e12t(:,:) * tms(:,:) ) 
     94      zbg_vfx_dyn = ztmp * glob_sum( wfx_dyn(:,:) * e12t(:,:) * tms(:,:) ) 
     95      zbg_vfx_bom = ztmp * glob_sum( wfx_bom(:,:) * e12t(:,:) * tms(:,:) ) 
     96      zbg_vfx_sum = ztmp * glob_sum( wfx_sum(:,:) * e12t(:,:) * tms(:,:) ) 
     97      zbg_vfx_res = ztmp * glob_sum( wfx_res(:,:) * e12t(:,:) * tms(:,:) ) 
     98      zbg_vfx_spr = ztmp * glob_sum( wfx_spr(:,:) * e12t(:,:) * tms(:,:) ) 
     99      zbg_vfx_snw = ztmp * glob_sum( wfx_snw(:,:) * e12t(:,:) * tms(:,:) ) 
     100      zbg_vfx_sub = ztmp * glob_sum( wfx_sub(:,:) * e12t(:,:) * tms(:,:) ) 
    101101 
    102102      ! Salt 
    103       zbg_sfx     = ztmp * glob_sum(     sfx(:,:) * area(:,:) * tms(:,:) ) 
    104       zbg_sfx_bri = ztmp * glob_sum( sfx_bri(:,:) * area(:,:) * tms(:,:) ) 
    105       zbg_sfx_res = ztmp * glob_sum( sfx_res(:,:) * area(:,:) * tms(:,:) ) 
    106       zbg_sfx_dyn = ztmp * glob_sum( sfx_dyn(:,:) * area(:,:) * tms(:,:) ) 
    107  
    108       zbg_sfx_bog = ztmp * glob_sum( sfx_bog(:,:) * area(:,:) * tms(:,:) ) 
    109       zbg_sfx_opw = ztmp * glob_sum( sfx_opw(:,:) * area(:,:) * tms(:,:) ) 
    110       zbg_sfx_sni = ztmp * glob_sum( sfx_sni(:,:) * area(:,:) * tms(:,:) ) 
    111       zbg_sfx_bom = ztmp * glob_sum( sfx_bom(:,:) * area(:,:) * tms(:,:) ) 
    112       zbg_sfx_sum = ztmp * glob_sum( sfx_sum(:,:) * area(:,:) * tms(:,:) ) 
     103      zbg_sfx     = ztmp * glob_sum(     sfx(:,:) * e12t(:,:) * tms(:,:) ) 
     104      zbg_sfx_bri = ztmp * glob_sum( sfx_bri(:,:) * e12t(:,:) * tms(:,:) ) 
     105      zbg_sfx_res = ztmp * glob_sum( sfx_res(:,:) * e12t(:,:) * tms(:,:) ) 
     106      zbg_sfx_dyn = ztmp * glob_sum( sfx_dyn(:,:) * e12t(:,:) * tms(:,:) ) 
     107 
     108      zbg_sfx_bog = ztmp * glob_sum( sfx_bog(:,:) * e12t(:,:) * tms(:,:) ) 
     109      zbg_sfx_opw = ztmp * glob_sum( sfx_opw(:,:) * e12t(:,:) * tms(:,:) ) 
     110      zbg_sfx_sni = ztmp * glob_sum( sfx_sni(:,:) * e12t(:,:) * tms(:,:) ) 
     111      zbg_sfx_bom = ztmp * glob_sum( sfx_bom(:,:) * e12t(:,:) * tms(:,:) ) 
     112      zbg_sfx_sum = ztmp * glob_sum( sfx_sum(:,:) * e12t(:,:) * tms(:,:) ) 
    113113 
    114114      ! Heat budget 
    115115      zbg_ihc      = glob_sum( et_i(:,:) * 1.e-20 ) ! ice heat content  [1.e-20 J] 
    116116      zbg_shc      = glob_sum( et_s(:,:) * 1.e-20 ) ! snow heat content [1.e-20 J] 
    117       zbg_hfx_dhc  = glob_sum( diag_heat_dhc(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
    118       zbg_hfx_spr  = glob_sum( hfx_spr(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
    119  
    120       zbg_hfx_thd  = glob_sum( hfx_thd(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
    121       zbg_hfx_dyn  = glob_sum( hfx_dyn(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
    122       zbg_hfx_res  = glob_sum( hfx_res(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
    123       zbg_hfx_sub  = glob_sum( hfx_sub(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
    124       zbg_hfx_snw  = glob_sum( hfx_snw(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
    125       zbg_hfx_sum  = glob_sum( hfx_sum(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
    126       zbg_hfx_bom  = glob_sum( hfx_bom(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
    127       zbg_hfx_bog  = glob_sum( hfx_bog(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
    128       zbg_hfx_dif  = glob_sum( hfx_dif(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
    129       zbg_hfx_opw  = glob_sum( hfx_opw(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
    130       zbg_hfx_out  = glob_sum( hfx_out(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
    131       zbg_hfx_in   = glob_sum(  hfx_in(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
     117      zbg_hfx_dhc  = glob_sum( diag_heat_dhc(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 
     118      zbg_hfx_spr  = glob_sum( hfx_spr(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 
     119 
     120      zbg_hfx_thd  = glob_sum( hfx_thd(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 
     121      zbg_hfx_dyn  = glob_sum( hfx_dyn(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 
     122      zbg_hfx_res  = glob_sum( hfx_res(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 
     123      zbg_hfx_sub  = glob_sum( hfx_sub(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 
     124      zbg_hfx_snw  = glob_sum( hfx_snw(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 
     125      zbg_hfx_sum  = glob_sum( hfx_sum(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 
     126      zbg_hfx_bom  = glob_sum( hfx_bom(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 
     127      zbg_hfx_bog  = glob_sum( hfx_bog(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 
     128      zbg_hfx_dif  = glob_sum( hfx_dif(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 
     129      zbg_hfx_opw  = glob_sum( hfx_opw(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 
     130      zbg_hfx_out  = glob_sum( hfx_out(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 
     131      zbg_hfx_in   = glob_sum(  hfx_in(:,:) * e12t(:,:) * tms(:,:) ) ! [in W] 
    132132     
    133133      ! --------------------------------------------- ! 
    134134      ! 2 - Trends due to forcing and ice growth/melt ! 
    135135      ! --------------------------------------------- ! 
    136       z_frc_vol = r1_rau0 * glob_sum( - emp(:,:) * area(:,:) * tms(:,:) ) ! volume fluxes 
    137       z_frc_sal = r1_rau0 * glob_sum(   sfx(:,:) * area(:,:) * tms(:,:) ) ! salt fluxes 
     136      z_frc_vol = r1_rau0 * glob_sum( - emp(:,:) * e12t(:,:) * tms(:,:) ) ! volume fluxes 
     137      z_frc_sal = r1_rau0 * glob_sum(   sfx(:,:) * e12t(:,:) * tms(:,:) ) ! salt fluxes 
    138138      z_bg_grme = glob_sum( - ( wfx_bog(:,:) + wfx_opw(:,:) + wfx_sni(:,:) + wfx_dyn(:,:) + & 
    139                           &     wfx_bom(:,:) + wfx_sum(:,:) + wfx_res(:,:) + wfx_snw(:,:) + wfx_sub(:,:) ) * area(:,:) * tms(:,:) ) ! volume fluxes 
     139                          &     wfx_bom(:,:) + wfx_sum(:,:) + wfx_res(:,:) + wfx_snw(:,:) + wfx_sub(:,:) ) * e12t(:,:) * tms(:,:) ) ! volume fluxes 
    140140      ! 
    141141      frc_vol  = frc_vol  + z_frc_vol  * rdt_ice 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90

    r5059 r5064  
    191191         CALL prt_ctl(tab2d_1=delta_i   , clinfo1=' lim_dyn  : delta_i   :') 
    192192         CALL prt_ctl(tab2d_1=strength  , clinfo1=' lim_dyn  : strength  :') 
    193          CALL prt_ctl(tab2d_1=area      , clinfo1=' lim_dyn  : cell area :') 
     193         CALL prt_ctl(tab2d_1=e12t      , clinfo1=' lim_dyn  : cell area :') 
    194194         CALL prt_ctl(tab2d_1=at_i      , clinfo1=' lim_dyn  : at_i      :') 
    195195         CALL prt_ctl(tab2d_1=vt_i      , clinfo1=' lim_dyn  : vt_i      :') 
     
    274274      rhoco  = rau0  * cw 
    275275      ! 
    276       !  Diffusion coefficients. 
     276      !  Diffusion coefficients 
    277277      zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 
    278278      IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain 
    279279            
    280280      za00 = ahi0 / ( 1.e05_wp )              ! 1.e05 = 100km = max grid space at 60° latitude in orca2 
    281                                               !                 (i.e. 60° = min latitude for ice cover)   
     281                                              !                    (60° = min latitude for ice cover)   
    282282      DO jj = 1, jpj 
    283283         DO ji = 1, jpi 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90

    r5055 r5064  
    345345                   ! Snow energy of melting 
    346346                   e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) 
    347                    ! Change dimensions 
    348                    e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac 
    349                    ! Multiply by volume, so that heat content in Joules 
    350                    e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * v_s(ji,jj,jl) / nlay_s 
     347 
     348                   ! Mutliply by volume, and divide by number of layers to get heat content in J/m2 
     349                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) / nlay_s 
    351350               END DO ! ji 
    352351            END DO ! jj 
     
    368367                      -   rcp     * ( ztmelts - rtt ) ) 
    369368 
    370                    ! Correct dimensions to avoid big values 
    371                    e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac  
    372  
    373                    ! Mutliply by ice volume, and divide by number of layers to get heat content in J 
    374                    e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / nlay_i 
     369                   ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 
     370                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) / nlay_i 
    375371               END DO ! ji 
    376372            END DO ! jj 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90

    r5055 r5064  
    339339            !-----------------------------------------------------------------------------! 
    340340            wfx_snw(ji,jj) = wfx_snw(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice     ! fresh water source for ocean 
    341             hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + esnow_mlt(ji,jj) * unit_fac / area(ji,jj) * r1_rdtice  ! heat sink for ocean (<0, W.m-2) 
     341            hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice     ! heat sink for ocean (<0, W.m-2) 
    342342 
    343343         END DO 
     
    381381         CALL prt_ctl_info(' - Cell values : ') 
    382382         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ') 
    383          CALL prt_ctl(tab2d_1=area , clinfo1=' lim_itd_me  : cell area :') 
     383         CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_itd_me  : cell area :') 
    384384         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_me  : at_i      :') 
    385385         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_me  : vt_i      :') 
     
    10931093               &                                + rhosn*vsrft(ji,jj)*(1.0-fsnowrft) 
    10941094 
    1095             ! in 1e-9 Joules (same as e_s) 
     1095            ! in J/m2 (same as e_s) 
    10961096            esnow_mlt(ji,jj) = esnow_mlt(ji,jj) - esrdg(ji,jj)*(1.0-fsnowrdg)         &   !rafting included 
    10971097               &                                - esrft(ji,jj)*(1.0-fsnowrft)           
     
    11331133               hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ji,jj,jk) * r1_rdtice  ! > 0 [W.m-2] ocean->ice flux  
    11341134 
    1135                ! Correct dimensions to avoid big values 
    1136                ersw(ji,jj,jk)   = ersw(ji,jj,jk) / unit_fac 
    1137  
    1138                ! Mutliply by ice volume, and divide by number of layers to get heat content in 1.e9 J 
    1139                ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean  
    1140                !! MV HC 2014 
    1141                ersw (ji,jj,jk)  = ersw(ji,jj,jk) * area(ji,jj) 
    1142  
     1135               ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean 
    11431136               erdg2(ji,jj,jk)  = erdg1(ji,jj,jk) + ersw(ji,jj,jk) 
    11441137 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90

    r5055 r5064  
    129129         DO jj = 1, jpj 
    130130            DO ji = 1, jpi 
    131                rswitch             = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) + epsi10 ) )     !0 if no ice and 1 if yes 
     131               rswitch           = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) + epsi10 ) )     !0 if no ice and 1 if yes 
    132132               ht_i(ji,jj,jl)    = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * rswitch 
    133                rswitch             = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i_b(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes 
     133               rswitch           = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i_b(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes 
    134134               zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * rswitch 
    135135               IF( a_i(ji,jj,jl) > epsi10 )   zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl)  
     
    167167      !----------------------------------------------------------------------------------------------- 
    168168      !- 4.1 Compute category boundaries 
    169       ! Tricky trick see limitd_me.F90 
    170       ! will be soon removed, CT 
    171       ! hi_max(kubnd) = 99. 
    172169      zhbnew(:,:,:) = 0._wp 
    173170 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limmsh.F90

    r5051 r5064  
    114114      CALL lbc_lnk( tmf(:,:), 'F', 1. )           ! lateral boundary conditions 
    115115 
    116       !                           !==  unmasked and masked area of T-grid cell 
    117       area(:,:) = e1t(:,:) * e2t(:,:) 
    118116      ! 
    119117   END SUBROUTINE lim_msh 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

    r5053 r5064  
    116116      CHARACTER (len=50) ::   charout 
    117117      REAL(wp) ::   zt11, zt12, zt21, zt22, ztagnx, ztagny, delta                         ! 
    118       REAL(wp) ::   za, zstms, zmask   ! local scalars 
    119       REAL(wp) ::   zc1, zc2, zc3             ! ice mass 
    120  
    121       REAL(wp) ::   dtevp              ! time step for subcycling 
    122       REAL(wp) ::   dtotel, ecc2, ecci ! square of yield ellipse eccenticity 
    123       REAL(wp) ::   z0, zr, zcca, zccb ! temporary scalars 
    124       REAL(wp) ::   zu_ice2, zv_ice1   ! 
    125       REAL(wp) ::   zddc, zdtc         ! delta on corners and on centre 
    126       REAL(wp) ::   zdst               ! shear at the center of the grid point 
    127       REAL(wp) ::   zdsshx, zdsshy     ! term for the gradient of ocean surface 
    128       REAL(wp) ::   sigma1, sigma2     ! internal ice stress 
     118      REAL(wp) ::   za, zstms          ! local scalars 
     119      REAL(wp) ::   zc1, zc2, zc3      ! ice mass 
     120 
     121      REAL(wp) ::   dtevp , z1_dtevp              ! time step for subcycling 
     122      REAL(wp) ::   dtotel, z1_dtotel, ecc2, ecci ! square of yield ellipse eccenticity 
     123      REAL(wp) ::   z0, zr, zcca, zccb            ! temporary scalars 
     124      REAL(wp) ::   zu_ice2, zv_ice1              ! 
     125      REAL(wp) ::   zddc, zdtc                    ! delta on corners and on centre 
     126      REAL(wp) ::   zdst                          ! shear at the center of the grid point 
     127      REAL(wp) ::   zdsshx, zdsshy                ! term for the gradient of ocean surface 
     128      REAL(wp) ::   sigma1, sigma2                ! internal ice stress 
    129129 
    130130      REAL(wp) ::   zresm         ! Maximal error on ice velocity 
     
    137137      REAL(wp), POINTER, DIMENSION(:,:) ::   zcorl1, zcorl2   ! coriolis parameter on U/V points 
    138138      REAL(wp), POINTER, DIMENSION(:,:) ::   za1ct , za2ct    ! temporary arrays 
    139       REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce1, v_oce1   ! ocean u/v component on U points                            
    140       REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce2, v_oce2   ! ocean u/v component on V points 
     139      REAL(wp), POINTER, DIMENSION(:,:) ::   v_oce1           ! ocean u/v component on U points                            
     140      REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce2           ! ocean u/v component on V points 
    141141      REAL(wp), POINTER, DIMENSION(:,:) ::   u_ice2, v_ice1   ! ice u/v component on V/U point 
    142142      REAL(wp), POINTER, DIMENSION(:,:) ::   zf1   , zf2      ! arrays for internal stresses 
     
    156156 
    157157      CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 
    158       CALL wrk_alloc( jpi,jpj, u_oce1, u_oce2, u_ice2, v_oce1 , v_oce2, v_ice1                ) 
     158      CALL wrk_alloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1                ) 
    159159      CALL wrk_alloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  ) 
    160160      CALL wrk_alloc( jpi,jpj, zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 ) 
     
    228228      !  zcorl2: Coriolis parameter on V-points                             
    229229      !  (ztagnx,ztagny): wind stress on U/V points                        
    230       !  u_oce1: ocean u component on u points                            
    231230      !  v_oce1: ocean v component on u points                           
    232231      !  u_oce2: ocean u component on v points                          
    233       !  v_oce2: ocean v component on v points                         
    234232 
    235233      IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==! 
     
    266264 
    267265            ! Mass, coriolis coeff. and currents 
    268             zmass1(ji,jj) = ( zt12*zc1 + zt11*zc2 ) / (zt11+zt12+zepsi) 
    269             zmass2(ji,jj) = ( zt22*zc1 + zt21*zc3 ) / (zt21+zt22+zepsi) 
    270             zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj)*fcor(ji,jj) + e1t(ji,jj)*fcor(ji+1,jj) )   & 
     266            zmass1(ji,jj) = ( zt12 * zc1 + zt11 * zc2 ) / ( zt11 + zt12 + zepsi ) 
     267            zmass2(ji,jj) = ( zt22 * zc1 + zt21 * zc3 ) / ( zt21 + zt22 + zepsi ) 
     268            zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj) * fcor(ji,jj) + e1t(ji,jj) * fcor(ji+1,jj) )   & 
    271269               &                          / ( e1t(ji,jj) + e1t(ji+1,jj) + zepsi ) 
    272             zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1)*fcor(ji,jj) + e2t(ji,jj)*fcor(ji,jj+1) )   & 
     270            zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1) * fcor(ji,jj) + e2t(ji,jj) * fcor(ji,jj+1) )   & 
    273271               &                          / ( e2t(ji,jj+1) + e2t(ji,jj) + zepsi ) 
    274272            ! 
    275             u_oce1(ji,jj)  = u_oce(ji,jj) 
    276             v_oce2(ji,jj)  = v_oce(ji,jj) 
    277  
    278273            ! Ocean has no slip boundary condition 
    279             v_oce1(ji,jj)  = 0.5*( (v_oce(ji,jj)+v_oce(ji,jj-1))*e1t(ji,jj)    & 
    280                &                 +(v_oce(ji+1,jj)+v_oce(ji+1,jj-1))*e1t(ji+1,jj)) & 
    281                &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)   
    282  
    283             u_oce2(ji,jj)  = 0.5*((u_oce(ji,jj)+u_oce(ji-1,jj))*e2t(ji,jj)     & 
    284                &                 +(u_oce(ji,jj+1)+u_oce(ji-1,jj+1))*e2t(ji,jj+1)) & 
    285                &                / (e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 
     274            v_oce1(ji,jj)  = 0.5 * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji,jj)      & 
     275               &                   + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji+1,jj) ) & 
     276               &                   / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj)   
     277 
     278            u_oce2(ji,jj)  = 0.5 * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj)      & 
     279               &                   + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj+1) ) & 
     280               &                   / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj) 
    286281 
    287282            ! Wind stress at U,V-point 
     
    316311      dtotel = dtevp / ( 2._wp * telast ) 
    317312#endif 
     313      z1_dtotel = 1._wp / ( 1._wp + dtotel ) 
     314      z1_dtevp  = 1._wp / dtevp 
    318315      !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter) 
    319316      ecc2 = ecc * ecc 
     
    334331 
    335332         DO jj = k_j1+1, k_jpj-1 
    336             DO ji = fs_2, jpim1   !RB bug no vect opt due to tmi 
     333            DO ji = fs_2, fs_jpim1   !RB bug no vect opt due to tmi 
    337334 
    338335               !   
     
    357354               divu_i(ji,jj) = (  e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
    358355                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
    359                   &            ) / area(ji,jj) 
     356                  &            ) * r1_e12t(ji,jj) 
    360357 
    361358               zdt(ji,jj) = ( ( u_ice(ji,jj) / e2u(ji,jj) - u_ice(ji-1,jj) / e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
    362359                  &         - ( v_ice(ji,jj) / e1v(ji,jj) - v_ice(ji,jj-1) / e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
    363                   &         ) / area(ji,jj) 
     360                  &         ) * r1_e12t(ji,jj) 
    364361 
    365362               ! 
    366363               zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) - u_ice(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   & 
    367364                  &         + ( v_ice(ji+1,jj) / e2v(ji+1,jj) - v_ice(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
    368                   &         ) / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2._wp - tmf(ji,jj) )   & 
     365                  &         ) * r1_e12f(ji,jj) * ( 2._wp - tmf(ji,jj) )   & 
    369366                  &         * tmi(ji,jj) * tmi(ji,jj+1) * tmi(ji+1,jj) * tmi(ji+1,jj+1) 
    370367 
    371368 
    372                v_ice1(ji,jj)  = 0.5_wp * ( ( v_ice(ji,jj)   + v_ice(ji,jj-1)   ) * e1t(ji+1,jj)   & 
    373                   &                      + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji,jj) )   & 
     369               v_ice1(ji,jj)  = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     & 
     370                  &                      + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) )   & 
    374371                  &                      / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj)  
    375372 
    376                u_ice2(ji,jj)  = 0.5_wp * ( ( u_ice(ji,jj)   + u_ice(ji-1,jj)   ) * e2t(ji,jj+1)   & 
    377                   &                      + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj) )   & 
     373               u_ice2(ji,jj)  = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
     374                  &                      + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) )   & 
    378375                  &                      / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj) 
    379  
    380             END DO 
    381          END DO 
    382          CALL lbc_lnk( v_ice1, 'U', -1. )   ;   CALL lbc_lnk( u_ice2, 'V', -1. )      ! lateral boundary cond. 
     376            END DO 
     377         END DO 
     378         CALL lbc_lnk( v_ice1  , 'U', -1. )   ;   CALL lbc_lnk( u_ice2  , 'V', -1. )      ! lateral boundary cond. 
    383379 
    384380         DO jj = k_j1+1, k_jpj-1 
     
    386382 
    387383               !- Calculate Delta at centre of grid cells 
    388                zdst      = ( e2u(ji, jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )   & 
    389                   &        + e1v(ji, jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1)   & 
    390                   &        ) / area(ji,jj) 
    391  
    392                delta = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst ) * usecc2 )   
     384               zdst          = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )   & 
     385                  &            + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1)   & 
     386                  &            ) * r1_e12t(ji,jj) 
     387 
     388               delta          = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 )   
    393389               delta_i(ji,jj) = delta + creepl 
    394                !-Calculate stress tensor components zs1 and zs2  
    395                !-at centre of grid cells (see section 3.5 of CICE user's guide). 
    396                zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( divu_i(ji,jj) / delta_i(ji,jj) - delta / delta_i(ji,jj) )  & 
    397                   &         * zpresh(ji,jj) ) ) / ( 1._wp + dtotel ) 
    398                zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / delta_i(ji,jj) * zpresh(ji,jj) ) )  & 
    399                   &         / ( 1._wp + dtotel ) 
    400  
    401             END DO 
    402          END DO 
    403  
    404          CALL lbc_lnk( zs1(:,:), 'T', 1. ) 
    405          CALL lbc_lnk( zs2(:,:), 'T', 1. ) 
    406  
    407          DO jj = k_j1+1, k_jpj-1 
    408             DO ji = fs_2, fs_jpim1 
     390 
    409391               !- Calculate Delta on corners 
    410                zddc  =      ( ( v_ice1(ji,jj+1)/e1u(ji,jj+1)                & 
    411                   &            -v_ice1(ji,jj)/e1u(ji,jj)                    & 
    412                   &           )*e1f(ji,jj)*e1f(ji,jj)                       & 
    413                   &          +( u_ice2(ji+1,jj)/e2v(ji+1,jj)                & 
    414                   &            -u_ice2(ji,jj)/e2v(ji,jj)                    & 
    415                   &           )*e2f(ji,jj)*e2f(ji,jj)                       & 
    416                   &         )                                               & 
    417                   &        / ( e1f(ji,jj) * e2f(ji,jj) ) 
    418  
    419                zdtc  =      (-( v_ice1(ji,jj+1)/e1u(ji,jj+1)                & 
    420                   &            -v_ice1(ji,jj)/e1u(ji,jj)                    & 
    421                   &           )*e1f(ji,jj)*e1f(ji,jj)                       & 
    422                   &          +( u_ice2(ji+1,jj)/e2v(ji+1,jj)                & 
    423                   &            -u_ice2(ji,jj)/e2v(ji,jj)                    & 
    424                   &           )*e2f(ji,jj)*e2f(ji,jj)                       & 
    425                   &         )                                               & 
    426                   &        / ( e1f(ji,jj) * e2f(ji,jj) ) 
    427  
    428                zddc = SQRT(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl 
    429  
    430                !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide). 
    431                zs12(ji,jj) = ( zs12(ji,jj) + dtotel *  & 
    432                   &          ( ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj) ) )  & 
    433                   &          / ( 1.0 + dtotel )  
    434  
    435             END DO ! ji 
    436          END DO ! jj 
    437  
    438          CALL lbc_lnk( zs12(:,:), 'F', 1. ) 
     392               zddc  = (  ( v_ice1(ji,jj+1) / e1u(ji,jj+1) - v_ice1(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  & 
     393                  &     + ( u_ice2(ji+1,jj) / e2v(ji+1,jj) - u_ice2(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  & 
     394                  &    ) * r1_e12f(ji,jj) 
     395 
     396               zdtc  = (- ( v_ice1(ji,jj+1) / e1u(ji,jj+1) - v_ice1(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  & 
     397                  &     + ( u_ice2(ji+1,jj) / e2v(ji+1,jj) - u_ice2(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  & 
     398                  &    ) * r1_e12f(ji,jj) 
     399 
     400               zddc = SQRT( zddc**2 + ( zdtc**2 + zds(ji,jj)**2 ) * usecc2 ) + creepl 
     401 
     402               !-Calculate stress tensor components zs1 and zs2 at centre of grid cells (see section 3.5 of CICE user's guide). 
     403               zs1(ji,jj)  = ( zs1 (ji,jj) + dtotel * ( divu_i(ji,jj) - delta ) / delta_i(ji,jj)   * zpresh(ji,jj)   & 
     404                  &          ) * z1_dtotel 
     405               zs2(ji,jj)  = ( zs2 (ji,jj) + dtotel *         ecci * zdt(ji,jj) / delta_i(ji,jj)   * zpresh(ji,jj)   & 
     406                  &          ) * z1_dtotel 
     407               !-Calculate stress tensor component zs12 at corners 
     408               zs12(ji,jj) = ( zs12(ji,jj) + dtotel *         ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj)  & 
     409                  &          ) * z1_dtotel  
     410 
     411            END DO 
     412         END DO 
     413         CALL lbc_lnk( zs1 , 'T', 1. )   ;   CALL lbc_lnk( zs2, 'T', 1. ) 
     414         CALL lbc_lnk( zs12, 'F', 1. ) 
    439415 
    440416         ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) 
     
    442418            DO ji = fs_2, fs_jpim1 
    443419               !- contribution of zs1, zs2 and zs12 to zf1 
    444                zf1(ji,jj) = 0.5*( (zs1(ji+1,jj)-zs1(ji,jj))*e2u(ji,jj) & 
    445                   &              +(zs2(ji+1,jj)*e2t(ji+1,jj)**2-zs2(ji,jj)*e2t(ji,jj)**2)/e2u(ji,jj) & 
    446                   &              +2.0*(zs12(ji,jj)*e1f(ji,jj)**2-zs12(ji,jj-1)*e1f(ji,jj-1)**2)/e1u(ji,jj) & 
    447                   &             ) / ( e1u(ji,jj)*e2u(ji,jj) ) 
     420               zf1(ji,jj) = 0.5 * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) & 
     421                  &             + ( zs2(ji+1,jj) * e2t(ji+1,jj)**2 - zs2(ji,jj) * e2t(ji,jj)**2 ) / e2u(ji,jj)          & 
     422                  &             + 2.0 * ( zs12(ji,jj) * e1f(ji,jj)**2 - zs12(ji,jj-1) * e1f(ji,jj-1)**2 ) / e1u(ji,jj) & 
     423                  &                ) * r1_e12u(ji,jj) 
    448424               ! contribution of zs1, zs2 and zs12 to zf2 
    449                zf2(ji,jj) = 0.5*( (zs1(ji,jj+1)-zs1(ji,jj))*e1v(ji,jj) & 
    450                   &              -(zs2(ji,jj+1)*e1t(ji,jj+1)**2 - zs2(ji,jj)*e1t(ji,jj)**2)/e1v(ji,jj) & 
    451                   &              + 2.0*(zs12(ji,jj)*e2f(ji,jj)**2 -    & 
    452                   zs12(ji-1,jj)*e2f(ji-1,jj)**2)/e2v(ji,jj) & 
    453                   &             ) / ( e1v(ji,jj)*e2v(ji,jj) ) 
     425               zf2(ji,jj) = 0.5 * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)  & 
     426                  &             - ( zs2(ji,jj+1) * e1t(ji,jj+1)**2 - zs2(ji,jj) * e1t(ji,jj)**2 ) / e1v(ji,jj)          & 
     427                  &             + 2.0 * ( zs12(ji,jj) * e2f(ji,jj)**2 - zs12(ji-1,jj) * e2f(ji-1,jj)**2 ) / e2v(ji,jj)  & 
     428                  &               )  * r1_e12v(ji,jj) 
    454429            END DO 
    455430         END DO 
     
    463438            DO jj = k_j1+1, k_jpj-1 
    464439               DO ji = fs_2, fs_jpim1 
    465                   zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj) 
    466                   z0           = zmass1(ji,jj)/dtevp 
     440                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * tmu(ji,jj) 
     441                  z0           = zmass1(ji,jj) * z1_dtevp 
    467442 
    468443                  ! SB modif because ocean has no slip boundary condition 
    469                   zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj)         & 
    470                      &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj))   & 
    471                      &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 
    472                   za           = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + & 
    473                      (zv_ice1-v_oce1(ji,jj))**2) * (1.0-zfrld1(ji,jj)) 
    474                   zr           = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & 
    475                      za*(u_oce1(ji,jj)) 
    476                   zcca         = z0+za 
     444                  zv_ice1      = 0.5 * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji  ,jj)     & 
     445                     &                 + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) )   & 
     446                     &                 / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj) 
     447                  za           = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 +  & 
     448                     &                         ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) ) 
     449                  zr           = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj) 
     450                  zcca         = z0 + za 
    477451                  zccb         = zcorl1(ji,jj) 
    478                   u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+zepsi)*zmask  
    479  
     452                  u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch  
    480453               END DO 
    481454            END DO 
     
    492465               DO ji = fs_2, fs_jpim1 
    493466 
    494                   zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj) 
    495                   z0           = zmass2(ji,jj)/dtevp 
     467                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * tmv(ji,jj) 
     468                  z0           = zmass2(ji,jj) * z1_dtevp 
    496469                  ! SB modif because ocean has no slip boundary condition 
    497                   zu_ice2       = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj)     & 
    498                      &                 + (u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1))   & 
    499                      &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 
    500                   za           = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + &  
    501                      (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 
    502                   zr           = z0*v_ice(ji,jj) + zf2(ji,jj) + & 
    503                      za2ct(ji,jj) + za*(v_oce2(ji,jj)) 
    504                   zcca         = z0+za 
     470                  zu_ice2      = 0.5 * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj)       & 
     471                     &                 + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) )   & 
     472                     &                 / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj) 
     473                  za           = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 +  &  
     474                     &                         ( v_ice(ji,jj) - v_oce(ji,jj))**2 ) * ( 1.0 - zfrld2(ji,jj) ) 
     475                  zr           = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj) 
     476                  zcca         = z0 + za 
    505477                  zccb         = zcorl2(ji,jj) 
    506                   v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+zepsi)*zmask 
    507  
     478                  v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch 
    508479               END DO 
    509480            END DO 
     
    520491            DO jj = k_j1+1, k_jpj-1 
    521492               DO ji = fs_2, fs_jpim1 
    522                   zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj) 
    523                   z0           = zmass2(ji,jj)/dtevp 
     493                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * tmv(ji,jj) 
     494                  z0           = zmass2(ji,jj) * z1_dtevp 
    524495                  ! SB modif because ocean has no slip boundary condition 
    525                   zu_ice2       = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj)      & 
    526                      &                 +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1))   & 
    527                      &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)    
    528  
    529                   za           = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + & 
    530                      (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 
    531                   zr           = z0*v_ice(ji,jj) + zf2(ji,jj) + & 
    532                      za2ct(ji,jj) + za*(v_oce2(ji,jj)) 
    533                   zcca         = z0+za 
     496                  zu_ice2      = 0.5 * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj)       & 
     497                     &                  +( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) )   & 
     498                     &                 / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj)    
     499 
     500                  za           = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 +  & 
     501                     &                         ( v_ice(ji,jj) - v_oce(ji,jj) )**2 ) * ( 1.0 - zfrld2(ji,jj) ) 
     502                  zr           = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj) 
     503                  zcca         = z0 + za 
    534504                  zccb         = zcorl2(ji,jj) 
    535                   v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+zepsi)*zmask 
    536  
     505                  v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch 
    537506               END DO 
    538507            END DO 
     
    548517            DO jj = k_j1+1, k_jpj-1 
    549518               DO ji = fs_2, fs_jpim1 
    550                   zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj) 
    551                   z0           = zmass1(ji,jj)/dtevp 
    552                   zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj)      & 
    553                      &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj))   & 
    554                      &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 
    555  
    556                   za           = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + & 
    557                      (zv_ice1-v_oce1(ji,jj))**2)*(1.0-zfrld1(ji,jj)) 
    558                   zr           = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & 
    559                      za*(u_oce1(ji,jj)) 
    560                   zcca         = z0+za 
     519                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * tmu(ji,jj) 
     520                  z0           = zmass1(ji,jj) * z1_dtevp 
     521                  zv_ice1      = 0.5 * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji,jj)       & 
     522                     &                 + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) )   & 
     523                     &                 / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj) 
     524 
     525                  za           = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 +  & 
     526                     &                         ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) ) 
     527                  zr           = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj) 
     528                  zcca         = z0 + za 
    561529                  zccb         = zcorl1(ji,jj) 
    562                   u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+zepsi)*zmask  
    563                END DO ! ji 
    564             END DO ! jj 
     530                  u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch  
     531               END DO 
     532            END DO 
    565533 
    566534            CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 
     
    577545            !---  Convergence test. 
    578546            DO jj = k_j1+1 , k_jpj-1 
    579                zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) ,           & 
    580                   ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 
    581             END DO 
    582             zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) ) 
     547               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 
     548            END DO 
     549            zresm = MAXVAL( zresr( 1:jpi, k_j1+1:k_jpj-1 ) ) 
    583550            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain 
    584551         ENDIF 
     
    637604            IF ( vt_i(ji,jj) <= zvmin ) THEN 
    638605 
    639                divu_i(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj)                      & 
    640                   &             -e2u(ji-1,jj)*u_ice(ji-1,jj)                  & 
    641                   &             +e1v(ji,jj)*v_ice(ji,jj)                      & 
    642                   &             -e1v(ji,jj-1)*v_ice(ji,jj-1)                  & 
    643                   &            )                                              & 
    644                   &            / area(ji,jj) 
    645  
    646                zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj)                    & 
    647                   &            -u_ice(ji-1,jj)/e2u(ji-1,jj)                & 
    648                   &           )*e2t(ji,jj)*e2t(ji,jj)                      & 
    649                   &          -( v_ice(ji,jj)/e1v(ji,jj)                    & 
    650                   &            -v_ice(ji,jj-1)/e1v(ji,jj-1)                & 
    651                   &           )*e1t(ji,jj)*e1t(ji,jj)                      & 
    652                   &         )                                              & 
    653                   &        / area(ji,jj) 
     606               divu_i(ji,jj) = (  e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj  ) * u_ice(ji-1,jj  )   & 
     607                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji  ,jj-1) * v_ice(ji  ,jj-1)   & 
     608                  &            ) * r1_e12t(ji,jj) 
     609 
     610               zdt(ji,jj) = ( ( u_ice(ji,jj) / e2u(ji,jj) - u_ice(ji-1,jj) / e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)  & 
     611                  &          -( v_ice(ji,jj) / e1v(ji,jj) - v_ice(ji,jj-1) / e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)  & 
     612                  &         ) * r1_e12t(ji,jj) 
    654613               ! 
    655614               ! SB modif because ocean has no slip boundary condition  
    656                zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1)              & 
    657                   &           - u_ice(ji,jj)   / e1u(ji,jj) )              & 
    658                   &           * e1f(ji,jj) * e1f(ji,jj)                    & 
    659                   &          + ( v_ice(ji+1,jj) / e2v(ji+1,jj)             & 
    660                   &            - v_ice(ji,jj)  / e2v(ji,jj) )              & 
    661                   &           * e2f(ji,jj) * e2f(ji,jj) )                  & 
    662                   &        / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) & 
    663                   &        * tmi(ji,jj) * tmi(ji,jj+1)                     & 
    664                   &        * tmi(ji+1,jj) * tmi(ji+1,jj+1) 
    665  
    666                zdst = (  e2u( ji  , jj   ) * v_ice1(ji  ,jj  )    & 
    667                   &           - e2u( ji-1, jj   ) * v_ice1(ji-1,jj  )    & 
    668                   &           + e1v( ji  , jj   ) * u_ice2(ji  ,jj  )    & 
    669                   &           - e1v( ji  , jj-1 ) * u_ice2(ji  ,jj-1)  ) / area(ji,jj) 
    670  
    671                delta = SQRT( divu_i(ji,jj)*divu_i(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 )   
     615               zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) - u_ice(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  & 
     616                  &          +( v_ice(ji+1,jj) / e2v(ji+1,jj) - v_ice(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  & 
     617                  &         ) * r1_e12f(ji,jj) * ( 2.0 - tmf(ji,jj) )                                     & 
     618                  &         * tmi(ji,jj) * tmi(ji,jj+1) * tmi(ji+1,jj) * tmi(ji+1,jj+1) 
     619 
     620               zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )    & 
     621                  &   + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1) ) * r1_e12t(ji,jj) 
     622 
     623               delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 )   
    672624               delta_i(ji,jj) = delta + creepl 
    673625             
    674626            ENDIF 
    675  
    676          END DO !jj 
    677       END DO !ji 
     627         END DO 
     628      END DO 
    678629      ! 
    679630      !------------------------------------------------------------------------------! 
     
    684635      DO jj = k_j1+1, k_jpj-1 
    685636         DO ji = fs_2, fs_jpim1 
    686             zdst= (  e2u( ji  , jj   ) * v_ice1(ji,jj)           &    
    687                &          - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)         &    
    688                &          + e1v( ji  , jj   ) * u_ice2(ji,jj)           &    
    689                &          - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1) ) / area(ji,jj)  
     637            zdst           = (  e2u(ji,jj) * v_ice1(ji,jj) - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)  &    
     638               &              + e1v(ji,jj) * u_ice2(ji,jj) - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1) ) * r1_e12t(ji,jj)  
    690639            shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst ) 
    691640         END DO 
     
    741690      ! 
    742691      CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 
    743       CALL wrk_dealloc( jpi,jpj, u_oce1, u_oce2, u_ice2, v_oce1 , v_oce2, v_ice1                ) 
     692      CALL wrk_dealloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1                ) 
    744693      CALL wrk_dealloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  ) 
    745694      CALL wrk_dealloc( jpi,jpj, zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 ) 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90

    r5055 r5064  
    2626   USE phycst           ! physical constants 
    2727   USE dom_oce          ! ocean domain 
    28    USE dom_ice,    ONLY : tms, area 
     28   USE dom_ice,    ONLY : tms 
    2929   USE ice              ! LIM sea-ice variables 
    3030   USE sbc_ice          ! Surface boundary condition: sea-ice fields 
     
    9999      !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108. 
    100100      !!              These refs are now obsolete since everything has been revised 
    101       !!              The ref should be Rousset et al., 2015? 
     101      !!              The ref should be Rousset et al., 2015 
    102102      !!--------------------------------------------------------------------- 
    103103      INTEGER, INTENT(in) ::   kt                                   ! number of iteration 
    104       ! 
    105104      INTEGER  ::   ji, jj, jl, jk                                  ! dummy loop indices 
    106       ! 
    107       REAL(wp) ::   zemp                                            !  local scalars 
     105      REAL(wp) ::   zemp                                            ! local scalars 
    108106      REAL(wp) ::   zf_mass                                         ! Heat flux associated with mass exchange ice->ocean (W.m-2) 
    109107      REAL(wp) ::   zfcm1                                           ! New solar flux received by the ocean 
     
    321319      !! ** input   : Namelist namicedia 
    322320      !!------------------------------------------------------------------- 
    323       REAL(wp) :: zsum, zarea 
    324       ! 
    325321      INTEGER  ::   ji, jj, jk                       ! dummy loop indices 
    326322      REAL(wp) ::   zcoefu, zcoefv, zcoeff          ! local scalar 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r5055 r5064  
    109109      ! 1.2) Heat content     
    110110      !-------------------- 
    111       ! Change the units of heat content; from global units to J.m3 
     111      ! Change the units of heat content; from J/m2 to J/m3 
    112112      DO jl = 1, jpl 
    113113         DO jk = 1, nlay_i 
     
    115115               DO ji = 1, jpi 
    116116                  !0 if no ice and 1 if yes 
    117                   rswitch = 1.0 - MAX(  0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 )  ) 
     117                  rswitch = 1.0 - MAX(  0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi20 )  ) 
    118118                  !Energy of melting q(S,T) [J.m-3] 
    119                   e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi10 ) ) * REAL( nlay_i ) 
    120                   !convert units ! very important that this line is here         
    121                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac  
     119                  e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL( nlay_i ) 
    122120               END DO 
    123121            END DO 
     
    127125               DO ji = 1, jpi 
    128126                  !0 if no ice and 1 if yes 
    129                   rswitch = 1.0 - MAX(  0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) + epsi10 )  ) 
     127                  rswitch = 1.0 - MAX(  0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) + epsi20 )  ) 
    130128                  !Energy of melting q(S,T) [J.m-3] 
    131                   e_s(ji,jj,jk,jl) = rswitch * e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL( nlay_s ) 
    132                   !convert units ! very important that this line is here 
    133                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * unit_fac  
     129                  e_s(ji,jj,jk,jl) = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL( nlay_s ) 
    134130               END DO 
    135131            END DO 
     
    453449      ! Ice heat content               
    454450      !------------------------ 
    455       ! Enthalpies are global variables we have to readjust the units (heat content in Joules) 
     451      ! Enthalpies are global variables we have to readjust the units (heat content in J/m2) 
    456452      DO jl = 1, jpl 
    457453         DO jk = 1, nlay_i 
    458             e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * a_i(:,:,jl) * ht_i(:,:,jl) / ( unit_fac * REAL( nlay_i ) ) 
     454            e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * a_i(:,:,jl) * ht_i(:,:,jl) / REAL( nlay_i ) 
    459455         END DO 
    460456      END DO 
     
    463459      ! Snow heat content               
    464460      !------------------------ 
    465       ! Enthalpies are global variables we have to readjust the units (heat content in Joules) 
     461      ! Enthalpies are global variables we have to readjust the units (heat content in J/m2) 
    466462      DO jl = 1, jpl 
    467463         DO jk = 1, nlay_s 
    468             e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * area(:,:) * a_i(:,:,jl) * ht_s(:,:,jl) / ( unit_fac * REAL( nlay_s ) ) 
     464            e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * a_i(:,:,jl) * ht_s(:,:,jl) / REAL( nlay_s ) 
    469465         END DO 
    470466      END DO 
     
    489485         CALL prt_ctl_info(' - Cell values : ') 
    490486         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ') 
    491          CALL prt_ctl(tab2d_1=area , clinfo1=' lim_thd  : cell area :') 
     487         CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_thd  : cell area :') 
    492488         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_thd  : at_i      :') 
    493489         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_thd  : vt_i      :') 
     
    520516      CALL wrk_dealloc( jpi, jpj, zqsr, zqns ) 
    521517 
     518      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    522519      !------------------------------------------------------------------------------| 
    523520      !  1) Transport of ice between thickness categories.                           | 
    524521      !------------------------------------------------------------------------------| 
    525       ! Given thermodynamic growth rates, transport ice between 
    526       ! thickness categories. 
     522      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     523 
     524      ! Given thermodynamic growth rates, transport ice between thickness categories. 
    527525      IF( jpl > 1 )   CALL lim_itd_th_rem( 1, jpl, kt ) 
    528526      ! 
     
    530528      CALL lim_var_agg(1) 
    531529 
     530      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     531      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    532532      !------------------------------------------------------------------------------| 
    533533      !  3) Add frazil ice growing in leads. 
     
    536536      CALL lim_var_glo2eqv    ! only for info 
    537537       
     538      ! conservation test 
     539      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     540 
    538541      IF(ln_ctl) THEN   ! Control print 
    539542         CALL prt_ctl_info(' ') 
    540543         CALL prt_ctl_info(' - Cell values : ') 
    541544         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ') 
    542          CALL prt_ctl(tab2d_1=area , clinfo1=' lim_itd_th  : cell area :') 
     545         CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_itd_th  : cell area :') 
    543546         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_th  : at_i      :') 
    544547         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_th  : vt_i      :') 
     
    568571      ENDIF 
    569572      ! 
    570       ! conservation test 
    571       IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    572  
    573573      IF( nn_timing == 1 )  CALL timing_stop('limthd') 
    574574 
     
    614614      !!                          ( dA = A/2h dh ) 
    615615      !!----------------------------------------------------------------------- 
    616       INTEGER, INTENT(in) ::   kideb, kiut   ! bounds for the spatial loop 
    617       INTEGER             ::   ji            ! dummy loop indices 
    618       REAL(wp)            ::   zhi_bef       ! ice thickness before thermo 
    619       REAL(wp)            ::   zdh_mel       ! net melting 
     616      INTEGER, INTENT(in) ::   kideb, kiut        ! bounds for the spatial loop 
     617      INTEGER             ::   ji                 ! dummy loop indices 
     618      REAL(wp)            ::   zhi_bef            ! ice thickness before thermo 
     619      REAL(wp)            ::   zdh_mel, zda_mel   ! net melting 
     620      REAL(wp)            ::   zv                 ! ice volume  
    620621 
    621622      DO ji = kideb, kiut 
    622          zhi_bef = ht_i_1d(ji) - ( dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) ) 
    623          zdh_mel = MIN( 0._wp,     dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) ) 
    624          IF( zdh_mel < 0._wp )   & 
    625             &   a_i_1d(ji) = MAX( 0._wp, a_i_1d(ji) + a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi10 ) ) )  
     623         zdh_mel = MIN( 0._wp, dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) ) 
     624         IF( zdh_mel < 0._wp )  THEN 
     625            zv         = a_i_1d(ji) * ht_i_1d(ji) 
     626            ! lateral melting = concentration change 
     627            zhi_bef     = ht_i_1d(ji) - zdh_mel 
     628            zda_mel     =  a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi10 ) ) 
     629            a_i_1d(ji)  = MAX( 0._wp, a_i_1d(ji) + zda_mel )  
     630            ! adjust thickness 
     631            rswitch     = 1._wp - MAX( 0._wp , SIGN( 1._wp , - a_i_1d(ji) + epsi20 ) ) 
     632            ht_i_1d(ji) = rswitch * zv / MAX( a_i_1d(ji), epsi20 ) 
     633            ! adjust thickness 
     634            at_i_1d(ji) = a_i_1d(ji) 
     635         END IF 
    626636      END DO 
    627       at_i_1d(:) = a_i_1d(:) 
    628637       
    629638   END SUBROUTINE lim_thd_lam 
     
    665674      IF ( ( jpl > 1 ) .AND. ( nn_monocat == 1 ) ) THEN  
    666675         nn_monocat = 0 
    667          WRITE(numout, *) ' nn_monocat must be 0 in multi-category case ' 
     676         IF(lwp) WRITE(numout, *) '  nn_monocat must be 0 in multi-category case ' 
    668677      ENDIF 
    669678 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90

    r5055 r5064  
    157157      REAL(wp)                            :: zh_thres      ! thickness thres. for G(h) computation 
    158158      REAL(wp)                            :: zhe           ! dummy factor 
    159       REAL(wp)                            :: zswitch       ! dummy switch 
    160159      REAL(wp)                            :: zkimean       ! mean sea ice thermal conductivity 
    161160      REAL(wp)                            :: zfac          ! dummy factor 
     
    363362         ! Fichefet and Morales Maqueda, JGR 1997 
    364363 
    365          zghe(:) = 0._wp 
     364         zghe(:) = 1._wp 
    366365 
    367366         SELECT CASE ( nn_monocat ) 
    368  
    369          CASE (0,2,4) 
    370  
    371             zghe(kideb:kiut) = 1._wp 
    372367 
    373368         CASE (1,3) ! LIM3 
     
    388383 
    389384               ! G(he) 
    390                zswitch  = MAX( 0._wp , SIGN( 1._wp , zhe - zh_thres ) )  ! =0 if zhe < zh_thres, if > 
    391                zghe(ji) = ( 1.0 - zswitch ) + zswitch * ( 0.5 + 0.5 * LOG( 2.*zhe / zepsilon ) ) 
     385               rswitch  = MAX( 0._wp , SIGN( 1._wp , zhe - zh_thres ) )  ! =0 if zhe < zh_thres, if > 
     386               zghe(ji) = ( 1._wp - rswitch ) + rswitch * 0.5_wp * ( 1._wp + LOG( 2.*zhe / zepsilon ) ) 
    392387    
    393388               ! Impose G(he) < 2. 
    394                zghe(ji) = MIN( zghe(ji), 2.0 ) 
     389               zghe(ji) = MIN( zghe(ji), 2._wp ) 
    395390 
    396391            END DO 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90

    r5055 r5064  
    132132      DO jk1 = 1, nlay_i 
    133133         DO ji = kideb, kiut 
    134             rswitch      = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zhnew(ji) + epsi10 ) )  
    135             qnew(ji,jk1) = rswitch * ( zqh_cum1(ji,jk1) - zqh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi10 ) 
     134            rswitch      = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zhnew(ji) + epsi20 ) )  
     135            qnew(ji,jk1) = rswitch * ( zqh_cum1(ji,jk1) - zqh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi20 ) 
    136136         ENDDO 
    137137      ENDDO 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90

    r5055 r5064  
    130130               DO ji = 1, jpi 
    131131                  !Energy of melting q(S,T) [J.m-3] 
    132                   rswitch          = 1._wp - MAX(  0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 )  )   !0 if no ice and 1 if yes 
    133                   e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) & 
    134                       &   / ( area(ji,jj) * MAX( v_i(ji,jj,jl) ,  epsi10 ) ) * REAL( nlay_i, wp ) 
    135                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac 
     132                  rswitch          = 1._wp - MAX(  0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi20 )  )   !0 if no ice 
     133                  e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl), epsi20 ) * REAL( nlay_i, wp ) 
    136134               END DO 
    137135            END DO 
     
    526524            DO jj = 1, jpj 
    527525               DO ji = 1, jpi 
    528                   ! heat content in Joules 
    529                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / ( REAL( nlay_i ,wp ) * unit_fac )  
     526                  ! heat content in J/m2 
     527                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) / REAL( nlay_i ,wp )  
    530528               END DO 
    531529            END DO 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90

    r5059 r5064  
    9595      ENDIF 
    9696 
    97       zsm(:,:) = area(:,:) 
     97      zsm(:,:) = e12t(:,:) 
    9898       
    9999      !                             !-------------------------------------! 
     
    153153         ! transported fields                                         
    154154         !------------------------- 
    155          zs0ow(:,:,1) = ato_i(:,:) * area(:,:)              ! Open water area  
    156          DO jl = 1, jpl 
    157             zs0sn (:,:,jl)   = v_s  (:,:,jl) * area(:,:)    ! Snow volume 
    158             zs0ice(:,:,jl)   = v_i  (:,:,jl) * area(:,:)    ! Ice  volume 
    159             zs0a  (:,:,jl)   = a_i  (:,:,jl) * area(:,:)    ! Ice area 
    160             zs0sm (:,:,jl)   = smv_i(:,:,jl) * area(:,:)    ! Salt content 
    161             zs0oi (:,:,jl)   = oa_i (:,:,jl) * area(:,:)    ! Age content 
    162             zs0c0 (:,:,jl)   = e_s  (:,:,1,jl)              ! Snow heat content 
    163             zs0e  (:,:,:,jl) = e_i  (:,:,:,jl)              ! Ice  heat content 
     155         zs0ow(:,:,1) = ato_i(:,:) * e12t(:,:)              ! Open water area  
     156         DO jl = 1, jpl 
     157            zs0sn (:,:,jl)   = v_s  (:,:,jl) * e12t(:,:)    ! Snow volume 
     158            zs0ice(:,:,jl)   = v_i  (:,:,jl) * e12t(:,:)    ! Ice  volume 
     159            zs0a  (:,:,jl)   = a_i  (:,:,jl) * e12t(:,:)    ! Ice area 
     160            zs0sm (:,:,jl)   = smv_i(:,:,jl) * e12t(:,:)    ! Salt content 
     161            zs0oi (:,:,jl)   = oa_i (:,:,jl) * e12t(:,:)    ! Age content 
     162            zs0c0 (:,:,jl)   = e_s  (:,:,1,jl) * e12t(:,:)  ! Snow heat content 
     163            DO jk = 1, nlay_i 
     164               zs0e  (:,:,jk,jl) = e_i  (:,:,jk,jl) * e12t(:,:) ! Ice  heat content 
     165            END DO 
    164166         END DO 
    165167 
     
    253255         ! Recover the properties from their contents 
    254256         !------------------------------------------- 
    255          ato_i(:,:) = zs0ow(:,:,1) / area(:,:) 
    256          DO jl = 1, jpl 
    257             v_i  (:,:,jl)   = zs0ice(:,:,jl) / area(:,:) 
    258             v_s  (:,:,jl)   = zs0sn (:,:,jl) / area(:,:) 
    259             smv_i(:,:,jl)   = zs0sm (:,:,jl) / area(:,:) 
    260             oa_i (:,:,jl)   = zs0oi (:,:,jl) / area(:,:) 
    261             a_i  (:,:,jl)   = zs0a  (:,:,jl) / area(:,:) 
    262             e_s  (:,:,1,jl) = zs0c0 (:,:,jl)               
    263             e_i  (:,:,:,jl) = zs0e  (:,:,:,jl)            
     257         ato_i(:,:) = zs0ow(:,:,1) * r1_e12t(:,:) 
     258         DO jl = 1, jpl 
     259            v_i  (:,:,jl)   = zs0ice(:,:,jl) * r1_e12t(:,:) 
     260            v_s  (:,:,jl)   = zs0sn (:,:,jl) * r1_e12t(:,:) 
     261            smv_i(:,:,jl)   = zs0sm (:,:,jl) * r1_e12t(:,:) 
     262            oa_i (:,:,jl)   = zs0oi (:,:,jl) * r1_e12t(:,:) 
     263            a_i  (:,:,jl)   = zs0a  (:,:,jl) * r1_e12t(:,:) 
     264            e_s  (:,:,1,jl) = zs0c0 (:,:,jl) * r1_e12t(:,:) 
     265            DO jk = 1, nlay_i 
     266               e_i  (:,:,jk,jl) = zs0e  (:,:,jk,jl) * r1_e12t(:,:) 
     267            END DO 
    264268         END DO 
    265269 
     
    380384                     wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 
    381385                     sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice  
    382                      hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 
    383                      hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 
     386                     hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0 
     387                     hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * r1_rdtice ! W.m-2 <0 
    384388                  ENDIF 
    385389 
     
    404408         DO jj = 1, jpj 
    405409            DO ji = 1, jpi 
    406                diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) / area(ji,jj) * unit_fac * r1_rdtice 
    407                diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) / area(ji,jj) * unit_fac * r1_rdtice 
     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 
    408412 
    409413               diag_trp_vi(ji,jj) = SUM( v_i(ji,jj,:) - zviold(ji,jj,:) ) * r1_rdtice 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90

    r5059 r5064  
    1616   USE sbc_ice         ! Surface boundary condition: ice fields 
    1717   USE dom_ice 
     18   USE dom_oce 
    1819   USE phycst          ! physical constants 
    1920   USE ice 
     
    148149            diag_heat_dhc(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  &  
    149150               &                     SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    & 
    150                &                   ) * unit_fac * r1_rdtice / area(ji,jj)    
     151               &                   ) * r1_rdtice    
    151152         END DO 
    152153      END DO 
     
    159160         CALL prt_ctl_info(' - Cell values : ') 
    160161         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ') 
    161          CALL prt_ctl(tab2d_1=area       , clinfo1=' lim_update1  : cell area   :') 
     162         CALL prt_ctl(tab2d_1=e12t       , clinfo1=' lim_update1  : cell area   :') 
    162163         CALL prt_ctl(tab2d_1=at_i       , clinfo1=' lim_update1  : at_i        :') 
    163164         CALL prt_ctl(tab2d_1=vt_i       , clinfo1=' lim_update1  : vt_i        :') 
     
    183184            CALL prt_ctl(tab2d_1=v_s        (:,:,jl)        , clinfo1= ' lim_update1  : v_s         : ') 
    184185            CALL prt_ctl(tab2d_1=v_s_b      (:,:,jl)        , clinfo1= ' lim_update1  : v_s_b       : ') 
    185             CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : e_i1        : ') 
    186             CALL prt_ctl(tab2d_1=e_i_b      (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : e_i1_b      : ') 
    187             CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : e_i2        : ') 
    188             CALL prt_ctl(tab2d_1=e_i_b      (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : e_i2_b      : ') 
     186            CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)      , clinfo1= ' lim_update1  : e_i1        : ') 
     187            CALL prt_ctl(tab2d_1=e_i_b      (:,:,1,jl)      , clinfo1= ' lim_update1  : e_i1_b      : ') 
     188            CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)      , clinfo1= ' lim_update1  : e_i2        : ') 
     189            CALL prt_ctl(tab2d_1=e_i_b      (:,:,2,jl)      , clinfo1= ' lim_update1  : e_i2_b      : ') 
    189190            CALL prt_ctl(tab2d_1=e_s        (:,:,1,jl)      , clinfo1= ' lim_update1  : e_snow      : ') 
    190191            CALL prt_ctl(tab2d_1=e_s_b      (:,:,1,jl)      , clinfo1= ' lim_update1  : e_snow_b    : ') 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90

    r5059 r5064  
    1616   USE sbc_ice         ! Surface boundary condition: ice fields 
    1717   USE dom_ice 
     18   USE dom_oce 
    1819   USE phycst          ! physical constants 
    1920   USE ice 
     
    190191               &                   ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  &  
    191192               &                     SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    & 
    192                &                   ) * unit_fac * r1_rdtice / area(ji,jj)    
     193               &                   ) * r1_rdtice    
    193194         END DO 
    194195      END DO 
     
    203204         CALL prt_ctl_info(' - Cell values : ') 
    204205         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ') 
    205          CALL prt_ctl(tab2d_1=area       , clinfo1=' lim_update2  : cell area   :') 
     206         CALL prt_ctl(tab2d_1=e12t       , clinfo1=' lim_update2  : cell area   :') 
    206207         CALL prt_ctl(tab2d_1=at_i       , clinfo1=' lim_update2  : at_i        :') 
    207208         CALL prt_ctl(tab2d_1=vt_i       , clinfo1=' lim_update2  : vt_i        :') 
     
    227228            CALL prt_ctl(tab2d_1=v_s        (:,:,jl)        , clinfo1= ' lim_update2  : v_s         : ') 
    228229            CALL prt_ctl(tab2d_1=v_s_b      (:,:,jl)        , clinfo1= ' lim_update2  : v_s_b       : ') 
    229             CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : e_i1        : ') 
    230             CALL prt_ctl(tab2d_1=e_i_b      (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : e_i1_b      : ') 
    231             CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2  : e_i2        : ') 
    232             CALL prt_ctl(tab2d_1=e_i_b      (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2  : e_i2_b      : ') 
     230            CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)      , clinfo1= ' lim_update2  : e_i1        : ') 
     231            CALL prt_ctl(tab2d_1=e_i_b      (:,:,1,jl)      , clinfo1= ' lim_update2  : e_i1_b      : ') 
     232            CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)      , clinfo1= ' lim_update2  : e_i2        : ') 
     233            CALL prt_ctl(tab2d_1=e_i_b      (:,:,2,jl)      , clinfo1= ' lim_update2  : e_i2_b      : ') 
    233234            CALL prt_ctl(tab2d_1=e_s        (:,:,1,jl)      , clinfo1= ' lim_update2  : e_snow      : ') 
    234235            CALL prt_ctl(tab2d_1=e_s_b      (:,:,1,jl)      , clinfo1= ' lim_update2  : e_snow_b    : ') 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90

    r5055 r5064  
    190190               DO ji = 1, jpi 
    191191                  !                                                              ! Energy of melting q(S,T) [J.m-3] 
    192                   rswitch   = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) )     ! rswitch = 0 if no ice and 1 if yes 
    193                   zq_i    = rswitch * e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi10 ) * REAL(nlay_i,wp)  
    194                   zq_i    = zq_i * unit_fac                             ! convert units 
     192                  rswitch   = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi20 ) )     ! rswitch = 0 if no ice and 1 if yes 
     193                  zq_i    = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL(nlay_i,wp)  
    195194                  ztmelts = -tmut * s_i(ji,jj,jk,jl) + rtt              ! Ice layer melt temperature 
    196195                  ! 
     
    216215               DO ji = 1, jpi 
    217216                  !Energy of melting q(S,T) [J.m-3] 
    218                   rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) + epsi10 ) )     ! rswitch = 0 if no ice and 1 if yes 
    219                   zq_s  = rswitch * e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL(nlay_s,wp) 
    220                   zq_s  = zq_s * unit_fac                                    ! convert units 
     217                  rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) + epsi20 ) )     ! rswitch = 0 if no ice and 1 if yes 
     218                  zq_s  = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL(nlay_s,wp) 
    221219                  ! 
    222220                  t_s(ji,jj,jk,jl) = rtt + rswitch * ( - zfac1 * zq_s + zfac2 ) 
     
    541539                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rtt * ( 1._wp - rswitch ) 
    542540                  ! update exchanges with ocean 
    543                   hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 
     541                  hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * r1_rdtice ! W.m-2 <0 
    544542               END DO 
    545543            END DO 
     
    579577               wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice 
    580578               wfx_snw(ji,jj)  = wfx_snw(ji,jj) - ( v_s(ji,jj,jl)   - zvs  ) * rhosn * r1_rdtice 
    581                hfx_res(ji,jj)  = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 
     579               hfx_res(ji,jj)  = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0 
    582580            END DO 
    583581         END DO 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90

    r5055 r5064  
    229229               ! Snow energy of melting 
    230230               e_s(ji,jj,jk,jl) = rswitch * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) 
    231                ! Change dimensions 
    232                e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac 
    233                ! Multiply by volume, so that heat content in 10^9 Joules 
    234                e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * v_s(ji,jj,jl) / nlay_s 
     231               ! Multiply by volume, so that heat content in J/m2 
     232               e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) / nlay_s 
    235233            END DO 
    236234            DO jk = 1, nlay_i 
     
    241239                  +   lfus    * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-epsi20) ) & 
    242240                  - rcp      * ( ztmelts - rtt ) ) 
    243                ! Correct dimensions to avoid big values 
    244                e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac  
    245                ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 
    246                e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / nlay_i 
     241               ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 
     242               e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / nlay_i 
    247243            END DO 
    248244 
Note: See TracChangeset for help on using the changeset viewer.