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 8239 for branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 – NEMO

Ignore:
Timestamp:
2017-06-28T17:55:50+02:00 (7 years ago)
Author:
clem
Message:

merge with v3_6_CMIP6_ice_diagnostics@r8238

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

    r7753 r8239  
    117117      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, zTauB, zTauE, zCor, zvel                          ! temporary scalars 
     119      REAL(wp) ::   zTauO, zTauB, zTauE, zvel                                ! temporary scalars 
    120120 
    121121      REAL(wp) ::   zsig1, zsig2                                             ! internal ice stress 
    122122      REAL(wp) ::   zresm                                                    ! Maximal error on ice velocity 
    123123      REAL(wp) ::   zintb, zintn                                             ! dummy argument 
     124      REAL(wp) ::   zfac_x, zfac_y 
    124125       
    125126      REAL(wp), POINTER, DIMENSION(:,:) ::   z1_e1t0, z1_e2t0                ! scale factors 
     
    140141                                                                             !   ocean surface (ssh_m) if ice is not embedded 
    141142                                                                             !   ice top surface if ice is embedded    
     143      REAL(wp), POINTER, DIMENSION(:,:) ::   zCorx, zCory                    ! Coriolis stress array 
     144      REAL(wp), POINTER, DIMENSION(:,:) ::   ztaux_oi, ztauy_oi              ! Ocean-to-ice stress array 
     145 
    142146      REAL(wp), POINTER, DIMENSION(:,:) ::   zswitchU, zswitchV              ! dummy arrays 
    143147      REAL(wp), POINTER, DIMENSION(:,:) ::   zmaskU, zmaskV                  ! mask for ice presence 
     
    154158      CALL wrk_alloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 
    155159      CALL wrk_alloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf ) 
     160      CALL wrk_alloc( jpi,jpj, zCorx, zCory) 
     161      CALL wrk_alloc( jpi,jpj, ztaux_oi, ztauy_oi) 
    156162 
    157163#if defined key_agrif  
     
    427433                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 
    428434 
     435                  ! Ocean-to-Ice stress 
     436                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
     437 
    429438                  ! tau_bottom/v_ice 
    430439                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 
     
    432441 
    433442                  ! Coriolis at V-points (energy conserving formulation) 
    434                   zCor  = - 0.25_wp * r1_e2v(ji,jj) *  & 
     443                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  & 
    435444                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
    436445                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
    437446 
    438447                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    439                   zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
     448                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
    440449 
    441450                  ! landfast switch => 0 = static friction ; 1 = sliding friction 
     
    451460                  ! Bouillon 2013 
    452461                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  & 
    453                   !!   &           + zfV(ji,jj) + zCor + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  & 
     462                  !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  & 
    454463                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj) 
    455464 
     
    471480                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
    472481 
     482                  ! Ocean-to-Ice stress 
     483                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
     484 
    473485                  ! tau_bottom/u_ice 
    474486                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) ) 
     
    476488 
    477489                  ! Coriolis at U-points (energy conserving formulation) 
    478                   zCor  =   0.25_wp * r1_e1u(ji,jj) *  & 
     490                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  & 
    479491                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
    480492                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
    481493                   
    482494                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    483                   zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
     495                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
    484496 
    485497                  ! landfast switch => 0 = static friction ; 1 = sliding friction 
     
    495507                  ! Bouillon 2013 
    496508                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  & 
    497                   !!   &           + zfU(ji,jj) + zCor + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  & 
     509                  !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  & 
    498510                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj) 
    499511               END DO 
     
    516528                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
    517529 
     530                  ! Ocean-to-Ice stress 
     531                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
     532 
    518533                  ! tau_bottom/u_ice 
    519534                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) ) 
     
    521536 
    522537                  ! Coriolis at U-points (energy conserving formulation) 
    523                   zCor  =   0.25_wp * r1_e1u(ji,jj) *  & 
     538                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  & 
    524539                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
    525540                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
    526541 
    527542                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    528                   zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
     543                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
    529544 
    530545                  ! landfast switch => 0 = static friction ; 1 = sliding friction 
     
    540555                  ! Bouillon 2013 
    541556                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  & 
    542                   !!   &           + zfU(ji,jj) + zCor + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  & 
     557                  !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  & 
    543558                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj) 
    544559               END DO 
     
    559574                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 
    560575 
     576                  ! Ocean-to-Ice stress 
     577                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
     578 
    561579                  ! tau_bottom/v_ice 
    562580                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 
     
    564582                   
    565583                  ! Coriolis at V-points (energy conserving formulation) 
    566                   zCor  = - 0.25_wp * r1_e2v(ji,jj) *  & 
     584                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  & 
    567585                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
    568586                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
    569587 
    570588                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    571                   zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
     589                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
    572590 
    573591                  ! landfast switch => 0 = static friction (tau_icebfr > zTauE); 1 = sliding friction 
     
    583601                  ! Bouillon 2013 
    584602                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  & 
    585                   !!   &           + zfV(ji,jj) + zCor + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  & 
     603                  !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  & 
    586604                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj) 
    587605               END DO 
     
    662680 
    663681      !------------------------------------------------------------------------------! 
    664       ! 5) Control prints of residual and charge ellipse 
     682      ! 5) SIMIP diagnostics 
     683      !------------------------------------------------------------------------------! 
     684                            
     685      DO jj = 2, jpjm1 
     686         DO ji = 2, jpim1 
     687             rswitch  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 
     688 
     689             ! Stress tensor invariants (normal and shear stress N/m) 
     690             diag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch                                 ! normal stress 
     691             diag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch   ! shear stress 
     692 
     693             ! Stress terms of the momentum equation (N/m2) 
     694             diag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch    ! sea surface slope stress term 
     695             diag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch 
     696 
     697             diag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch    ! Coriolis stress term 
     698             diag_corstry(ji,jj) = zCory(ji,jj) * rswitch 
     699 
     700             diag_intstrx(ji,jj) = zfU(ji,jj)   * rswitch    ! internal stress term 
     701             diag_intstry(ji,jj) = zfV(ji,jj)   * rswitch 
     702            
     703             diag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch  ! oceanic stress 
     704             diag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch 
     705 
     706             ! 2D ice mass, snow mass, area transport arrays (X, Y) 
     707             zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch 
     708             zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch 
     709 
     710             diag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component 
     711             diag_ymtrp_ice(ji,jj) = rhoic * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   '' 
     712 
     713             diag_xmtrp_snw(ji,jj) = rhosn * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component 
     714             diag_ymtrp_snw(ji,jj) = rhosn * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   '' 
     715 
     716             diag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )         ! area transport,      X-component 
     717             diag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )         !        ''            Y-   '' 
     718 
     719         END DO 
     720      END DO 
     721 
     722      CALL lbc_lnk_multi(   diag_sig1   , 'T',  1., diag_sig2   , 'T',  1.,   & 
     723                 &          diag_dssh_dx, 'U', -1., diag_dssh_dy, 'V', -1.,   & 
     724                 &          diag_corstrx, 'U', -1., diag_corstry, 'V', -1.,   &  
     725                 &          diag_intstrx, 'U', -1., diag_intstry, 'V', -1.    ) 
     726 
     727      CALL lbc_lnk_multi(   diag_utau_oi, 'U', -1., diag_vtau_oi, 'V', -1.    ) 
     728 
     729      ! 
     730      !------------------------------------------------------------------------------! 
     731      ! 6) Control prints of residual and charge ellipse 
    665732      !------------------------------------------------------------------------------! 
    666733      ! 
     
    703770      CALL wrk_dealloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 
    704771      CALL wrk_dealloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf ) 
     772      CALL wrk_dealloc( jpi,jpj, zCorx, zCory ) 
     773      CALL wrk_dealloc( jpi,jpj, ztaux_oi, ztauy_oi ) 
    705774 
    706775   END SUBROUTINE lim_rhg 
Note: See TracChangeset for help on using the changeset viewer.