Changeset 10652


Ignore:
Timestamp:
2019-02-08T11:29:16+01:00 (20 months ago)
Author:
mathiot
Message:

branch and suggested changes for ticket #2229 (iceberg melt along coastlines)

Location:
NEMO/branches/2019/fix_ticket2229
Files:
1 edited
1 copied

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/fix_ticket2229/src/OCE/ICB/icbutl.F90

    r10570 r10652  
    131131      !!             is half the off shore value, wile the normal-to-the-coast value is zero. 
    132132      !!             This is OK as a starting point. 
     133      !!       !!pm  HARD CODED: - rho_air now computed in sbcblk (what are the effect ?) 
     134      !!                         - drag coefficient (should it be namelist parameter ?) 
    133135      !! 
    134136      !!---------------------------------------------------------------------- 
     
    142144      !!---------------------------------------------------------------------- 
    143145 
    144       pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj )     ! scale factors 
     146      pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj )      ! scale factors 
    145147      pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj ) 
    146148      ! 
    147       puo  = icb_utl_bilin_h( uo_e, pi, pj, 'U' )             ! ocean velocities 
    148       pvo  = icb_utl_bilin_h( vo_e, pi, pj, 'V' ) 
    149       psst = icb_utl_bilin_h( tt_e, pi, pj, 'T' )             ! SST 
    150       pcn  = icb_utl_bilin_h( fr_e , pi, pj, 'T' )            ! ice concentration 
    151       pff  = icb_utl_bilin_h( ff_e , pi, pj, 'F' )            ! Coriolis parameter 
    152       ! 
    153       pua  = icb_utl_bilin_h( ua_e , pi, pj, 'U' )            ! 10m wind 
    154       pva  = icb_utl_bilin_h( va_e , pi, pj, 'V' )            ! here (ua,va) are stress => rough conversion from stress to speed 
    155       zcd  = 1.22_wp * 1.5e-3_wp                              ! air density * drag coefficient 
     149      puo  = icb_utl_bilin_h( uo_e, pi, pj, 'U', .false.  )    ! ocean velocities 
     150      pvo  = icb_utl_bilin_h( vo_e, pi, pj, 'V', .false. ) 
     151      psst = icb_utl_bilin_h( tt_e, pi, pj, 'T', .true.   )    ! SST 
     152      pcn  = icb_utl_bilin_h( fr_e, pi, pj, 'T', .true.   )    ! ice concentration 
     153      pff  = icb_utl_bilin_h( ff_e, pi, pj, 'F', .false.  )    ! Coriolis parameter 
     154      ! 
     155      pua  = icb_utl_bilin_h( ua_e, pi, pj, 'U', .true.   )    ! 10m wind 
     156      pva  = icb_utl_bilin_h( va_e, pi, pj, 'V', .true.   )    ! here (ua,va) are stress => rough conversion from stress to speed 
     157      zcd  = 1.22_wp * 1.5e-3_wp                               ! air density * drag coefficient  
    156158      zmod = 1._wp / MAX(  1.e-20, SQRT(  zcd * SQRT( pua*pua + pva*pva)  )  ) 
    157159      pua  = pua * zmod                                       ! note: stress module=0 necessarly implies ua=va=0 
     
    159161 
    160162#if defined key_si3 
    161       pui = icb_utl_bilin_h( ui_e , pi, pj, 'U' )              ! sea-ice velocities 
    162       pvi = icb_utl_bilin_h( vi_e , pi, pj, 'V' ) 
    163       phi = icb_utl_bilin_h( hicth, pi, pj, 'T' )              ! ice thickness 
     163      pui = icb_utl_bilin_h( ui_e , pi, pj, 'U', .false. )    ! sea-ice velocities 
     164      pvi = icb_utl_bilin_h( vi_e , pi, pj, 'V', .false. ) 
     165      phi = icb_utl_bilin_h( hicth, pi, pj, 'T', .true.  )    ! ice thickness 
    164166#else 
    165167      pui = 0._wp 
     
    169171 
    170172      ! Estimate SSH gradient in i- and j-direction (centred evaluation) 
    171       pssh_i = ( icb_utl_bilin_h( ssh_e, pi+0.1_wp, pj, 'T' ) -   & 
    172          &       icb_utl_bilin_h( ssh_e, pi-0.1_wp, pj, 'T' )  ) / ( 0.2_wp * pe1 ) 
    173       pssh_j = ( icb_utl_bilin_h( ssh_e, pi, pj+0.1_wp, 'T' ) -   & 
    174          &       icb_utl_bilin_h( ssh_e, pi, pj-0.1_wp, 'T' )  ) / ( 0.2_wp * pe2 ) 
     173      pssh_i = ( icb_utl_bilin_h( ssh_e, pi+0.1_wp, pj, 'T', .true. ) -   & 
     174         &       icb_utl_bilin_h( ssh_e, pi-0.1_wp, pj, 'T', .true. )  ) / ( 0.2_wp * pe1 ) 
     175      pssh_j = ( icb_utl_bilin_h( ssh_e, pi, pj+0.1_wp, 'T', .true. ) -   & 
     176         &       icb_utl_bilin_h( ssh_e, pi, pj-0.1_wp, 'T', .true. )  ) / ( 0.2_wp * pe2 ) 
    175177      ! 
    176178   END SUBROUTINE icb_utl_interp 
    177179 
    178180 
    179    REAL(wp) FUNCTION icb_utl_bilin_h( pfld, pi, pj, cd_type ) 
     181   REAL(wp) FUNCTION icb_utl_bilin_h( pfld, pi, pj, cd_type, plmask ) 
    180182      !!---------------------------------------------------------------------- 
    181183      !!                  ***  FUNCTION icb_utl_bilin  *** 
     
    191193      REAL(wp)                            , INTENT(in) ::   pi, pj    ! targeted coordinates in (i,j) referential 
    192194      CHARACTER(len=1)                    , INTENT(in) ::   cd_type   ! type of pfld array grid-points: = T , U , V or F points 
     195      LOGICAL                             , INTENT(in) ::   plmask    ! special treatment of mask point 
    193196      ! 
    194197      INTEGER  ::   ii, ij   ! local integer 
    195198      REAL(wp) ::   zi, zj   ! local real 
     199      REAL(wp) :: zw1, zw2, zw3, zw4 
     200      REAL(wp), DIMENSION(4) :: zmask 
    196201      !!---------------------------------------------------------------------- 
    197202      ! 
     
    224229      ! find position in this processor. Prevent near edge problems (see #1389) 
    225230      ! 
    226       IF    ( ii < mig( 1 ) ) THEN   ;   ii = 1 
    227       ELSEIF( ii > mig(jpi) ) THEN   ;   ii = jpi 
     231      IF    ( ii < mig( 1 ) ) THEN   ;   ii = 1   ;  
     232      ELSEIF( ii > mig(jpi) ) THEN   ;   ii = jpi ; 
    228233      ELSE                           ;   ii = mi1(ii) 
    229234      ENDIF 
    230       IF    ( ij < mjg( 1 ) ) THEN   ;   ij = 1 
    231       ELSEIF( ij > mjg(jpj) ) THEN   ;   ij = jpj 
     235      IF    ( ij < mjg( 1 ) ) THEN   ;   ij = 1   ; 
     236      ELSEIF( ij > mjg(jpj) ) THEN   ;   ij = jpj ; 
    232237      ELSE                           ;   ij  = mj1(ij) 
    233238      ENDIF 
    234239      ! 
    235       IF( ii == jpi )   ii = ii-1       
    236       IF( ij == jpj )   ij = ij-1 
    237       ! 
    238       icb_utl_bilin_h = ( pfld(ii,ij  ) * (1.-zi) + pfld(ii+1,ij  ) * zi ) * (1.-zj)   & 
    239          &            + ( pfld(ii,ij+1) * (1.-zi) + pfld(ii+1,ij+1) * zi ) *     zj 
     240      IF( ii == jpi ) ii = ii-1 
     241      IF( ij == jpj ) ij = ij-1 
     242      ! 
     243      ! define mask array  
     244      IF (plmask) THEN 
     245         ! land value is not used in the interpolation 
     246         SELECT CASE ( cd_type ) 
     247         CASE ( 'T' ) 
     248            zmask = (/tmask(ii,ij,1), tmask(ii+1,ij,1), tmask(ii,ij+1,1), tmask(ii+1,ij+1,1)/) 
     249         CASE ( 'U' ) 
     250            zmask = (/umask(ii,ij,1), umask(ii+1,ij,1), umask(ii,ij+1,1), umask(ii+1,ij+1,1)/) 
     251         CASE ( 'V' ) 
     252            zmask = (/vmask(ii,ij,1), vmask(ii+1,ij,1), vmask(ii,ij+1,1), vmask(ii+1,ij+1,1)/) 
     253         CASE ( 'F' ) 
     254            ! F case only used for coriolis, ff_f is not mask so zmask = 1 
     255            zmask = 1. 
     256         END SELECT 
     257      ELSE 
     258         ! land value is used during interpolation 
     259         zmask = 1. 
     260      END iF 
     261      ! 
     262      ! compute weight 
     263      zw1 = zmask(1) * (1._wp-zi) * (1._wp-zj) 
     264      zw2 = zmask(2) *        zi  * (1._wp-zj) 
     265      zw3 = zmask(3) * (1._wp-zi) *        zj 
     266      zw4 = zmask(4) *        zi  *        zj 
     267      ! 
     268      ! compute interpolated value 
     269      icb_utl_bilin_h = ( pfld(ii,ij)*zw1 + pfld(ii+1,ij)*zw2 + pfld(ii,ij+1)*zw3 + pfld(ii+1,ij+1)*zw4 ) / MAX(1.e-20, zw1+zw2+zw3+zw4)  
    240270      ! 
    241271   END FUNCTION icb_utl_bilin_h 
Note: See TracChangeset for help on using the changeset viewer.