Changeset 10691


Ignore:
Timestamp:
2019-02-15T17:41:15+01:00 (18 months ago)
Author:
mathiot
Message:

merge fix_ticket2229 into the trunk (ticket #2229)

Location:
NEMO/trunk/src/OCE/ICB
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/OCE/ICB/icb_oce.F90

    r10425 r10691  
    8989   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ua_e, va_e 
    9090   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ssh_e 
     91   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   tmask_e, umask_e, vmask_e 
    9192#if defined key_si3 || defined key_cice 
    9293   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ui_e, vi_e 
     
    169170      ! 
    170171      ! expanded arrays for bilinear interpolation 
    171       ALLOCATE( uo_e(0:jpi+1,0:jpj+1) , ua_e(0:jpi+1,0:jpj+1) ,   & 
    172          &      vo_e(0:jpi+1,0:jpj+1) , va_e(0:jpi+1,0:jpj+1) ,   & 
     172      ALLOCATE( uo_e(0:jpi+1,0:jpj+1) , ua_e(0:jpi+1,0:jpj+1) ,    & 
     173         &      vo_e(0:jpi+1,0:jpj+1) , va_e(0:jpi+1,0:jpj+1) ,    & 
    173174#if defined key_si3 || defined key_cice 
    174175         &      ui_e(0:jpi+1,0:jpj+1) ,                            & 
     
    183184      icb_alloc = icb_alloc + ill 
    184185 
     186      ALLOCATE( tmask_e(0:jpi+1,0:jpj+1), umask_e(0:jpi+1,0:jpj+1), vmask_e(0:jpi+1,0:jpj+1), & 
     187         &      STAT=ill) 
     188      icb_alloc = icb_alloc + ill 
     189 
    185190      ALLOCATE( nicbfldpts(jpi) , nicbflddest(jpi) , nicbfldproc(jpni) , & 
    186191         &      nicbfldnsend(jpni), nicbfldexpect(jpni) , nicbfldreq(jpni), STAT=ill) 
  • NEMO/trunk/src/OCE/ICB/icbini.F90

    r10570 r10691  
    224224      src_calving_hflx(:,:) = 0._wp 
    225225 
     226      ! definition of extended surface masked needed by icb_bilin_h 
     227      tmask_e(:,:) = 0._wp   ;   tmask_e(1:jpi,1:jpj) = tmask(:,:,1) 
     228      umask_e(:,:) = 0._wp   ;   umask_e(1:jpi,1:jpj) = umask(:,:,1) 
     229      vmask_e(:,:) = 0._wp   ;   vmask_e(1:jpi,1:jpj) = vmask(:,:,1) 
     230      CALL lbc_lnk_icb( 'icbini', tmask_e, 'T', +1._wp, 1, 1 ) 
     231      CALL lbc_lnk_icb( 'icbini', umask_e, 'T', +1._wp, 1, 1 ) 
     232      CALL lbc_lnk_icb( 'icbini', vmask_e, 'T', +1._wp, 1, 1 ) 
     233      ! 
    226234      ! assign each new iceberg with a unique number constructed from the processor number 
    227235      ! and incremented by the total number of processors 
  • NEMO/trunk/src/OCE/ICB/icbutl.F90

    r10570 r10691  
    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_e(ii,ij), tmask_e(ii+1,ij), tmask_e(ii,ij+1), tmask_e(ii+1,ij+1)/) 
     249         CASE ( 'U' ) 
     250            zmask = (/umask_e(ii,ij), umask_e(ii+1,ij), umask_e(ii,ij+1), umask_e(ii+1,ij+1)/) 
     251         CASE ( 'V' ) 
     252            zmask = (/vmask_e(ii,ij), vmask_e(ii+1,ij), vmask_e(ii,ij+1), vmask_e(ii+1,ij+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.