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

Changeset 6796


Ignore:
Timestamp:
2016-07-06T17:35:54+02:00 (8 years ago)
Author:
clem
Message:

rewriting of ice rheology to conserve energy

Location:
branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref

    r6774 r6796  
    9393                                    !                      a very large value ensures ice velocity=0 even with a small contact area 
    9494                                    !                      recommended range: ?? (should be greater than atm-ice stress => >0.1 N/m2) 
     95   rn_lfrelax     =    1.e-5        !     (ln_landfast=T)  relaxation time scale to reach static friction (s-1)                  
    9596/ 
    9697!------------------------------------------------------------------------------ 

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

    r6763 r6796  
    215215   REAL(wp), PUBLIC ::   rn_gamma         !: fraction of ocean depth that ice must reach to initiate landfast ice 
    216216   REAL(wp), PUBLIC ::   rn_icebfr        !: maximum bottom stress per unit area of contact (landfast ice)  
     217   REAL(wp), PUBLIC ::   rn_lfrelax       !: relaxation time scale (s-1) to reach static friction (landfast ice)  
    217218 
    218219   !                                     !!** ice-diffusion namelist (namicehdf) ** 

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

    r6746 r6796  
    7878      IF( ln_landfast ) THEN 
    7979         DO jl = 1, jpl 
    80             WHERE( ht_i(:,:,jl) > bathy(:,:) * rn_gamma )  tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 
     80            WHERE( ht_i(:,:,jl) > ht(:,:) * rn_gamma )  tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 
    8181         END DO 
    8282      ENDIF 
    … …  
    148148      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    149149      NAMELIST/namicedyn/ nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio,  rn_creepl, rn_ecc, & 
    150          &                nn_nevp, rn_relast, ln_landfast, rn_gamma, rn_icebfr 
     150         &                nn_nevp, rn_relast, ln_landfast, rn_gamma, rn_icebfr, rn_lfrelax 
    151151      !!------------------------------------------------------------------- 
    152152 
    … …  
    177177         WRITE(numout,*) '   Landfast: fraction of ocean depth that ice must reach       rn_gamma      = ', rn_gamma 
    178178         WRITE(numout,*) '   Landfast: maximum bottom stress per unit area of contact    rn_icebfr     = ', rn_icebfr 
     179         WRITE(numout,*) '   Landfast: relax time scale (s-1) to reach static friction   rn_lfrelax    = ', rn_lfrelax 
    179180      ENDIF 
    180181      ! 

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

    r6746 r6796  
    847847         END DO 
    848848    
    849          strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:) 
     849         strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:) * tmask(:,:,1) 
    850850                         ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation 
    851851         ksmooth = 1 
    … …  
    856856      ELSE                      ! kstrngth ne 1:  Hibler (1979) form 
    857857         ! 
    858          strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) )  ) 
     858         strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) )  ) * tmask(:,:,1) 
    859859         ! 
    860860         ksmooth = 1 

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

    r6789 r6796  
    1010   !!            3.4  !  2011-01  (A. Porter)  dynamical allocation  
    1111   !!            3.5  !  2012-08  (R. Benshila)  AGRIF 
    12    !!            3.6  !  2016-06  (C. Rousset) Rewriting and add the possibility to use mEVP from Bouillon 2013 
     12   !!            3.6  !  2016-06  (C. Rousset) Rewriting + landfast ice + possibility to use mEVP (Bouillon 2013) 
    1313   !!---------------------------------------------------------------------- 
    1414#if defined key_lim3 
    … …  
    112112      CHARACTER (len=50) ::   charout 
    113113 
    114       REAL(wp) ::   dtevp, z1_dtevp                                          ! time step for subcycling 
     114      REAL(wp) ::   zdtevp, z1_dtevp                                         ! time step for subcycling 
    115115      REAL(wp) ::   ecc2, z1_ecc2                                            ! square of yield ellipse eccenticity 
    116116      REAL(wp) ::   zbeta, zalph1, z1_alph1, zalph2, z1_alph2                ! alpha and beta from Bouillon 2009 and 2013 
    117       REAL(wp) ::   zm1, zm2, zm3                                            ! ice/snow mass 
     117      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                            ! ice/snow mass 
    118118      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2            ! temporary scalars 
    119       REAL(wp) ::   zTauO, zTauA, zTauB, ZCor, zmt, zu_ice2, zv_ice1, zvel   ! temporary scalars 
     119      REAL(wp) ::   zTauO, zTauB, zTauE, zCor, zvel                          ! temporary scalars 
    120120 
    121121      REAL(wp) ::   zsig1, zsig2                                             ! internal ice stress 
    … …  
    123123      REAL(wp) ::   zintb, zintn                                             ! dummy argument 
    124124       
    125       REAL(wp), POINTER, DIMENSION(:,:) ::   zpresh                          ! temporary array for ice strength 
    126125      REAL(wp), POINTER, DIMENSION(:,:) ::   z1_e1t0, z1_e2t0                ! scale factors 
    127126      REAL(wp), POINTER, DIMENSION(:,:) ::   zp_delt                         ! P/delta at T points 
    128127      ! 
    129       REAL(wp), POINTER, DIMENSION(:,:) ::   za1   , za2                     ! ice fraction on U/V points 
    130       REAL(wp), POINTER, DIMENSION(:,:) ::   zmass1, zmass2                  ! ice/snow mass on U/V points 
    131       REAL(wp), POINTER, DIMENSION(:,:) ::   zcor1 , zcor2                   ! coriolis parameter on U/V points 
    132       REAL(wp), POINTER, DIMENSION(:,:) ::   zspgu , zspgv                   ! surface pressure gradient at U/V points 
    133       REAL(wp), POINTER, DIMENSION(:,:) ::   v_oce1, u_oce2, v_ice1, u_ice2  ! ocean/ice u/v component on V/U points                            
    134       REAL(wp), POINTER, DIMENSION(:,:) ::   zf1   , zf2                     ! internal stresses 
     128      REAL(wp), POINTER, DIMENSION(:,:) ::   zaU   , zaV                     ! ice fraction on U/V points 
     129      REAL(wp), POINTER, DIMENSION(:,:) ::   zmU_t, zmV_t                    ! ice/snow mass/dt on U/V points 
     130      REAL(wp), POINTER, DIMENSION(:,:) ::   zmf                             ! coriolis parameter at T points 
     131      REAL(wp), POINTER, DIMENSION(:,:) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points 
     132      REAL(wp), POINTER, DIMENSION(:,:) ::   zspgU , zspgV                   ! surface pressure gradient at U/V points 
     133      REAL(wp), POINTER, DIMENSION(:,:) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                            
     134      REAL(wp), POINTER, DIMENSION(:,:) ::   zfU   , zfV                     ! internal stresses 
    135135       
    136136      REAL(wp), POINTER, DIMENSION(:,:) ::   zds                             ! shear 
    … …  
    148148      !!------------------------------------------------------------------- 
    149149 
    150       CALL wrk_alloc( jpi,jpj, zpresh, z1_e1t0, z1_e2t0, zp_delt ) 
    151       CALL wrk_alloc( jpi,jpj, za1, za2, zmass1, zmass2, zcor1, zcor2 ) 
    152       CALL wrk_alloc( jpi,jpj, zspgu, zspgv, v_oce1, u_oce2, v_ice1, u_ice2, zf1, zf2 ) 
     150      CALL wrk_alloc( jpi,jpj, z1_e1t0, z1_e2t0, zp_delt ) 
     151      CALL wrk_alloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia ) 
     152      CALL wrk_alloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV ) 
    153153      CALL wrk_alloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 
    154154      CALL wrk_alloc( jpi,jpj, zswitch1, zswitch2, zfmask, zwf ) 
    … …  
    205205 
    206206      ! Time step for subcycling 
    207       dtevp    = rdt_ice / REAL( nn_nevp ) 
    208       z1_dtevp = 1._wp / dtevp 
     207      zdtevp   = rdt_ice / REAL( nn_nevp ) 
     208      z1_dtevp = 1._wp / zdtevp 
    209209 
    210210      ! alpha parameters (Bouillon 2009) 
    … …  
    228228      ! Ice strength 
    229229      CALL lim_itd_me_icestrength( nn_icestr ) 
    230       zpresh(:,:) = tmask(:,:,1) *  strength(:,:) 
    231230 
    232231      ! scale factors 
    … …  
    263262 
    264263            ! ice fraction at U-V points 
    265             za1(ji,jj) = ( at_i(ji,jj) * e1t(ji+1,jj) + at_i(ji+1,jj) * e1t(ji,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
    266             za2(ji,jj) = ( at_i(ji,jj) * e2t(ji,jj+1) + at_i(ji,jj+1) * e2t(ji,jj) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
     264            zaU(ji,jj) = ( at_i(ji,jj) * e1t(ji+1,jj) + at_i(ji+1,jj) * e1t(ji,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
     265            zaV(ji,jj) = ( at_i(ji,jj) * e2t(ji,jj+1) + at_i(ji,jj+1) * e2t(ji,jj) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
    267266 
    268267            ! Ice/snow mass at U-V points 
    … …  
    270269            zm2 = ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) ) 
    271270            zm3 = ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) ) 
    272             zmass1(ji,jj) = ( zm1 * e1t(ji+1,jj) + zm2 * e1t(ji,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
    273             zmass2(ji,jj) = ( zm1 * e2t(ji,jj+1) + zm3 * e2t(ji,jj) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
     271            zmassU = ( zm1 * e1t(ji+1,jj) + zm2 * e1t(ji,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
     272            zmassV = ( zm1 * e2t(ji,jj+1) + zm3 * e2t(ji,jj) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
    274273 
    275274            ! Ocean currents at U-V points 
    276             v_oce1(ji,jj)  = 0.5_wp * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji+1,jj)    & 
    277                &                      + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
    278              
    279             u_oce2(ji,jj)  = 0.5_wp * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj+1)    & 
    280                &                      + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
    281  
    282             ! Coriolis at U-V points (m*f) 
    283             zcor1(ji,jj) =   zmass1(ji,jj) * 0.5_wp * ( ff(ji,jj) + ff(ji  ,jj-1) ) 
    284             zcor2(ji,jj) = - zmass2(ji,jj) * 0.5_wp * ( ff(ji,jj) + ff(ji-1,jj  ) ) 
     275            v_oceU(ji,jj)   = 0.5_wp * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji+1,jj)    & 
     276               &                       + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
     277             
     278            u_oceV(ji,jj)   = 0.5_wp * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj+1)    & 
     279               &                       + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
     280 
     281            ! Coriolis at T points (m*f) 
     282            zmf(ji,jj)      = zm1 * ff_t(ji,jj) 
     283 
     284            ! m/dt 
     285            zmU_t(ji,jj)    = zmassU * z1_dtevp 
     286            zmV_t(ji,jj)    = zmassV * z1_dtevp 
     287 
     288            ! Drag ice-atm. 
     289            zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) 
     290            zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj) 
    285291 
    286292            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points 
    287             zspgu(ji,jj) = - zmass1(ji,jj) * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj) 
    288             zspgv(ji,jj) = - zmass2(ji,jj) * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj) 
     293            zspgU(ji,jj)    = - zmassU * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj) 
     294            zspgV(ji,jj)    = - zmassV * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj) 
    289295 
    290296            ! switches 
    291             zswitch1(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) 
    292             zswitch2(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) 
     297            zswitch1(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) ) 
     298            zswitch2(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) ) 
    293299 
    294300         END DO 
    295301      END DO 
     302      CALL lbc_lnk( zmf, 'T', 1. ) 
    296303      ! 
    297304      !------------------------------------------------------------------------------! 
    … …  
    311318         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 
    312319         DO jj = 1, jpjm1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points 
    313             DO ji = 1, fs_jpim1 
     320            DO ji = 1, jpim1 
    314321 
    315322               ! shear at F points 
    … …  
    346353 
    347354               ! P/delta at T points 
    348                zp_delt(ji,jj) = zpresh(ji,jj) / ( zdelta + rn_creepl ) 
     355               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 
    349356                
    350357               ! stress at T points 
    … …  
    375382 
    376383               ! U points 
    377                zf1(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             & 
     384               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             & 
    378385                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    & 
    379386                  &                    ) * r1_e2u(ji,jj)                                                                      & 
    … …  
    383390 
    384391               ! V points 
    385                zf2(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             & 
     392               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             & 
    386393                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    & 
    387394                  &                    ) * r1_e1v(ji,jj)                                                                      & 
    … …  
    391398 
    392399               ! u_ice at V point 
    393                u_ice2(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
     400               u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
    394401                  &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
    395402                
    396403               ! v_ice at U point 
    397                v_ice1(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     & 
     404               v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     & 
    398405                  &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
    399406 
    … …  
    402409         ! 
    403410         ! --- Computation of ice velocity --- ! 
    404          !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen smart and large nn_nevp 
     411         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp 
    405412         !  Bouillon et al. 2009 (eq 34-35) => stable 
    406413         IF( MOD(jter,2) .EQ. 0 ) THEN ! even iterations 
    … …  
    409416               DO ji = fs_2, fs_jpim1 
    410417 
    411                   zvel    = SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_ice2(ji,jj) * u_ice2(ji,jj) ) 
     418                  ! tau_io/(v_oce - v_ice) 
     419                  zTauO = zaV(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
     420                     &                             + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 
     421 
     422                  ! tau_bottom/v_ice 
     423                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 
     424                  ztauB = - tau_icebfr(ji,jj) / zvel 
     425 
     426                  ! Coriolis at V-points (energy conserving formulation) 
     427                  zCor  = - 0.25_wp * r1_e2v(ji,jj) *  & 
     428                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
     429                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
     430 
     431                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     432                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
     433 
     434                  ! landfast switch => 0 = static friction ; 1 = sliding friction 
     435                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztauB * v_ice(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
    412436                   
    413                   zTauA   =   za2(ji,jj) * vtau_ice(ji,jj) 
    414                   zTauO   =   za2(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
    415                      &                                 + ( u_ice2(ji,jj) - u_oce2(ji,jj) ) * ( u_ice2(ji,jj) - u_oce2(ji,jj) ) ) 
    416                   ztauB   = - tau_icebfr(ji,jj) / MAX( zvel, zepsi ) 
    417                   zCor    =   zcor2(ji,jj)  * u_ice2(ji,jj) 
    418                   zmt     =   zmass2(ji,jj) * z1_dtevp 
    419                    
    420                   ! Bouillon 2009 
    421                   v_ice(ji,jj) = ( zmt * v_ice(ji,jj) + zf2(ji,jj) + zCor + zTauA + zTauO * v_oce(ji,jj) + zspgv(ji,jj)  & 
    422                      &           ) / MAX( zmt + zTauO - zTauB, zepsi ) * zswitch2(ji,jj) 
     437                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
     438                  v_ice(ji,jj) = (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity 
     439                     &                                   + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     440                     &                                   ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast 
     441                     &           + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
     442                     &           ) * zswitch2(ji,jj) 
    423443                  ! Bouillon 2013 
    424                   !!v_ice(ji,jj) = ( zmt * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  & 
    425                   !!   &           + zf2(ji,jj) + zCor + zTauA + zTauO * v_oce(ji,jj) + zspgv(ji,jj)  & 
    426                   !!   &           ) / MAX( zmt * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch2(ji,jj) 
     444                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  & 
     445                  !!   &           + zfV(ji,jj) + zCor + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  & 
     446                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch2(ji,jj) 
    427447 
    428448               END DO 
    … …  
    439459            DO jj = 2, jpjm1 
    440460               DO ji = fs_2, fs_jpim1 
    441  
    442                   ! v_ice at U point 
    443                   zv_ice1 = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     & 
    444                      &               + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)  
     461                                
     462                  ! tau_io/(u_oce - u_ice) 
     463                  zTauO = zaU(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
     464                     &                             + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
     465 
     466                  ! tau_bottom/u_ice 
     467                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) ) 
     468                  ztauB = - tau_icebfr(ji,jj) / zvel 
     469 
     470                  ! Coriolis at U-points (energy conserving formulation) 
     471                  zCor  =   0.25_wp * r1_e1u(ji,jj) *  & 
     472                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
     473                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
    445474                   
    446                   zvel    = SQRT( v_ice1(ji,jj) * v_ice1(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 
    447                    
    448                   zTauA   =   za1(ji,jj) * utau_ice(ji,jj) 
    449                   zTauO   =   za1(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
    450                      &                                 + ( v_ice1(ji,jj) - v_oce1(ji,jj) ) * ( v_ice1(ji,jj) - v_oce1(ji,jj) ) ) 
    451                   ztauB   = - tau_icebfr(ji,jj) / MAX( zvel, zepsi ) 
    452                   zCor    =   zcor1(ji,jj)  * zv_ice1 
    453                   zmt     =   zmass1(ji,jj) * z1_dtevp 
    454                    
    455                   ! Bouillon 2009 
    456                   u_ice(ji,jj) = ( zmt * u_ice(ji,jj) + zf1(ji,jj) + zCor + zTauA + zTauO * u_oce(ji,jj) + zspgu(ji,jj)  & 
    457                      &           ) / MAX( zmt + zTauO - zTauB, zepsi ) * zswitch1(ji,jj) 
     475                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     476                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
     477 
     478                  ! landfast switch => 0 = static friction ; 1 = sliding friction 
     479                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztauB * u_ice(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     480 
     481                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
     482                  u_ice(ji,jj) = (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity 
     483                     &                                   + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     484                     &                                   ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast 
     485                     &           + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
     486                     &           ) * zswitch1(ji,jj) 
    458487                  ! Bouillon 2013 
    459                   !!u_ice(ji,jj) = ( zmt * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  & 
    460                   !!   &           + zf1(ji,jj) + zCor + zTauA + zTauO * u_oce(ji,jj) + zspgu(ji,jj)  & 
    461                   !!   &           ) / MAX( zmt * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch1(ji,jj) 
     488                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  & 
     489                  !!   &           + zfU(ji,jj) + zCor + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  & 
     490                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch1(ji,jj) 
    462491               END DO 
    463492            END DO 
    … …  
    476505               DO ji = fs_2, fs_jpim1 
    477506 
    478                   zvel    = SQRT( v_ice1(ji,jj) * v_ice1(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 
    479  
    480                   zTauA   =   za1(ji,jj) * utau_ice(ji,jj) 
    481                   zTauO   =   za1(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
    482                      &                                 + ( v_ice1(ji,jj) - v_oce1(ji,jj) ) * ( v_ice1(ji,jj) - v_oce1(ji,jj) ) ) 
    483                   ztauB   = - tau_icebfr(ji,jj) / MAX( zvel, zepsi ) 
    484                   zCor    =   zcor1(ji,jj)  * v_ice1(ji,jj) 
    485                   zmt     =   zmass1(ji,jj) * z1_dtevp 
    486  
    487                   ! Bouillon 2009 
    488                   u_ice(ji,jj) = ( zmt * u_ice(ji,jj) + zf1(ji,jj) + zCor + zTauA + zTauO * u_oce(ji,jj) + zspgu(ji,jj)  & 
    489                      &           ) / MAX( zmt + zTauO - zTauB, zepsi ) * zswitch1(ji,jj) 
     507                  ! tau_io/(u_oce - u_ice) 
     508                  zTauO = zaU(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
     509                     &                             + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
     510 
     511                  ! tau_bottom/u_ice 
     512                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) ) 
     513                  ztauB = - tau_icebfr(ji,jj) / zvel 
     514 
     515                  ! Coriolis at U-points (energy conserving formulation) 
     516                  zCor  =   0.25_wp * r1_e1u(ji,jj) *  & 
     517                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
     518                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
     519 
     520                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     521                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
     522 
     523                  ! landfast switch => 0 = static friction ; 1 = sliding friction 
     524                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztauB * u_ice(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     525 
     526                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
     527                  u_ice(ji,jj) = (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity 
     528                     &                                   + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     529                     &                                   ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )            &  ! m/dt + tau_io(only ice part) + landfast 
     530                     &           + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
     531                     &           ) * zswitch1(ji,jj) 
    490532                  ! Bouillon 2013 
    491                   !!u_ice(ji,jj) = ( zmt * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  & 
    492                   !!   &           + zf1(ji,jj) + zCor + zTauA + zTauO * u_oce(ji,jj) + zspgu(ji,jj)  & 
    493                   !!   &           ) / MAX( zmt * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch1(ji,jj) 
     533                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  & 
     534                  !!   &           + zfU(ji,jj) + zCor + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  & 
     535                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch1(ji,jj) 
    494536               END DO 
    495537            END DO 
    … …  
    505547            DO jj = 2, jpjm1 
    506548               DO ji = fs_2, fs_jpim1 
    507  
    508                   ! u_ice at V point 
    509                   zu_ice2 = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
    510                      &               + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
    511  
    512                   zvel    = SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_ice2(ji,jj) * u_ice2(ji,jj) ) 
    513549          
    514                   zTauA   =   za2(ji,jj) * vtau_ice(ji,jj) 
    515                   zTauO   =   za2(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
    516                      &                                 + ( u_ice2(ji,jj) - u_oce2(ji,jj) ) * ( u_ice2(ji,jj) - u_oce2(ji,jj) ) ) 
    517                   ztauB   = - tau_icebfr(ji,jj) / MAX( zvel, zepsi ) 
    518                   zCor    =   zcor2(ji,jj)  * zu_ice2 
    519                   zmt     =   zmass2(ji,jj) * z1_dtevp 
     550                  ! tau_io/(v_oce - v_ice) 
     551                  zTauO = zaV(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
     552                     &                             + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 
     553 
     554                  ! tau_bottom/v_ice 
     555                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 
     556                  ztauB = - tau_icebfr(ji,jj) / zvel 
    520557                   
    521                   ! Bouillon 2009 
    522                   v_ice(ji,jj) = ( zmt * v_ice(ji,jj) + zf2(ji,jj) + zCor + zTauA + zTauO * v_oce(ji,jj) + zspgv(ji,jj)  & 
    523                      &           ) / MAX( zmt + zTauO - zTauB, zepsi ) * zswitch2(ji,jj) 
     558                  ! Coriolis at V-points (energy conserving formulation) 
     559                  zCor  = - 0.25_wp * r1_e2v(ji,jj) *  & 
     560                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
     561                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
     562 
     563                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     564                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
     565 
     566                  ! landfast switch => 0 = static friction ; 1 = sliding friction 
     567                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztauB * v_ice(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     568                   
     569                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
     570                  v_ice(ji,jj) = (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity 
     571                     &                                   + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     572                     &                                   ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )            &  ! m/dt + tau_io(only ice part) + landfast 
     573                     &           + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
     574                     &           ) * zswitch2(ji,jj) 
    524575                  ! Bouillon 2013 
    525                   !!v_ice(ji,jj) = ( zmt * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  & 
    526                   !!   &           + zf2(ji,jj) + zCor + zTauA + zTauO * v_oce(ji,jj) + zspgv(ji,jj)  & 
    527                   !!   &           ) / MAX( zmt * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch2(ji,jj) 
     576                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  & 
     577                  !!   &           + zfV(ji,jj) + zCor + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  & 
     578                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch2(ji,jj) 
    528579               END DO 
    529580            END DO 
    … …  
    579630      !------------------------------------------------------------------------------! 
    580631      DO jj = 1, jpjm1 
    581          DO ji = 1, fs_jpim1 
     632         DO ji = 1, jpim1 
    582633 
    583634            ! shear at F points 
    … …  
    649700            DO jj = 2, jpjm1 
    650701               DO ji = 2, jpim1 
    651                   IF (zpresh(ji,jj) > 1.0) THEN 
    652                      zsig1 = ( zs1(ji,jj) + SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*zpresh(ji,jj) )  
    653                      zsig2 = ( zs1(ji,jj) - SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*zpresh(ji,jj) ) 
     702                  IF (strength(ji,jj) > 1.0) THEN 
     703                     zsig1 = ( zs1(ji,jj) + SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*strength(ji,jj) )  
     704                     zsig2 = ( zs1(ji,jj) - SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*strength(ji,jj) ) 
    654705                     WRITE(charout,FMT="('lim_rhg  :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)") 
    655706                     CALL prt_ctl_info(charout) 
    … …  
    662713      ENDIF 
    663714      ! 
    664       CALL wrk_dealloc( jpi,jpj, zpresh, z1_e1t0, z1_e2t0, zp_delt ) 
    665       CALL wrk_dealloc( jpi,jpj, za1, za2, zmass1, zmass2, zcor1, zcor2 ) 
    666       CALL wrk_dealloc( jpi,jpj, zspgu, zspgv, v_oce1, u_oce2, v_ice1, u_ice2, zf1, zf2 ) 
     715      CALL wrk_dealloc( jpi,jpj, z1_e1t0, z1_e2t0, zp_delt ) 
     716      CALL wrk_dealloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia ) 
     717      CALL wrk_dealloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV ) 
    667718      CALL wrk_dealloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 
    668719      CALL wrk_dealloc( jpi,jpj, zswitch1, zswitch2, zfmask, zwf ) 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r5123 r6796  
    168168   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  e1e2t          !: surface at t-point (m2) 
    169169   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ff             !: coriolis factor (2.*omega*sin(yphi) ) (s-1) 
     170   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ff_t           !: clem: coriolis factor at T-point (s-1) 
    170171 
    171172   !!---------------------------------------------------------------------- 
    … …  
    350351         &      glamv(jpi,jpj) , gphiv(jpi,jpj) , e1v(jpi,jpj) , e2v(jpi,jpj) , r1_e1v(jpi,jpj) , r1_e2v(jpi,jpj) ,   &   
    351352         &      glamf(jpi,jpj) , gphif(jpi,jpj) , e1f(jpi,jpj) , e2f(jpi,jpj) , r1_e1f(jpi,jpj) , r1_e2f(jpi,jpj) ,   & 
    352          &      e1e2t(jpi,jpj) , ff  (jpi,jpj) , STAT=ierr(3) )      
     353         &      e1e2t(jpi,jpj) , ff_t (jpi,jpj) , ff (jpi,jpj) , STAT=ierr(3) )      
    353354         ! 
    354355      ALLOCATE( gdep3w_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0 (jpi,jpj,jpk) ,                         & 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90

    r6204 r6796  
    528528      CASE ( 0, 1, 4 )               ! mesh on the sphere 
    529529 
    530          ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) )  
     530         ff  (:,:) = 2. * omega * SIN( rad * gphif(:,:) )  
     531         ff_t(:,:) = 2. * omega * SIN( rad * gphit(:,:) )  ! clem: coriolis at T-point  
    531532 
    532533      CASE ( 2 )                     ! f-plane at ppgphi0  
    533534 
    534          ff(:,:) = 2. * omega * SIN( rad * ppgphi0 ) 
     535         ff  (:,:) = 2. * omega * SIN( rad * ppgphi0 ) 
     536         ff_t(:,:) = 2. * omega * SIN( rad * ppgphi0 )  ! clem: coriolis at T-point 
    535537 
    536538         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff(1,1) 
    … …  
    551553         zf0     = 2. * omega * SIN( rad * zphi0 )                              ! compute f0 1st point south 
    552554 
    553          ff(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south) 
     555         ff  (:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south) 
     556         ff_t(:,:) = ( zf0  + zbeta * gphit(:,:) * 1.e+3 )                        ! clem: coriolis at T-point 
    554557          
    555558         IF(lwp) THEN 
    … …  
    572575         zf0   = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south 
    573576 
    574          ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south) 
     577         ff  (:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south) 
     578         ff_t(:,:) = ( zf0 + zbeta * ABS( gphit(:,:) - zphi0 ) * rad * ra )   ! clem: coriolis at T-point 
    575579 
    576580         IF(lwp) THEN 

  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/TOOLS/REBUILD_NEMO/icb_combrest.py

    r6449 r6796  
    169169    sys.exit(15) 
    170170  fo = Dataset(pathout, 'w') 
    171   for dim in ['x','y','c']: 
     171  for dim in ['x','y','c','k']: 
    172172    indim = fi.dimensions[dim] 
    173173    fo.createDimension(dim, len(indim)) 
    174   for var in ['calving','calving_hflx','stored_ice','stored_heat']: 
     174  for var in ['kount','calving','calving_hflx','stored_ice','stored_heat']: 
    175175    invar = fi.variables[var] 
    176176    fo.createVariable(var, invar.datatype, invar.dimensions) 
    177177    fo.variables[var][:] = invar[:] 
    178     fo.variables[var].long_name = invar.long_name 
    179     fo.variables[var].units = invar.units 
     178    if "long_name" in invar.ncattrs(): 
     179        fo.variables[var].long_name = invar.long_name 
     180    if "units" in invar.ncattrs(): 
     181        fo.variables[var].units = invar.units 
    180182  os.remove(pathout.replace('.nc','_WORK.nc')) 
    181183# 
Note: See TracChangeset for help on using the changeset viewer.