Changeset 13265


Ignore:
Timestamp:
2020-07-08T13:34:28+02:00 (5 weeks ago)
Author:
mathiot
Message:

ticket #2494 and #2375: ticket #2494 changes and part of ticket #2375 (lbc_icb_lnk on extended variable at T point removed), ff_e initialisation and lbc_lnk move in the icbini.F90; tests on possible useless lbc_lnk in icbclv not yet done

Location:
NEMO/branches/2020/tickets_2494_2375/src/OCE/ICB
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/tickets_2494_2375/src/OCE/ICB/icb_oce.F90

    r12472 r13265  
    8686   ! particularly for MPP when iceberg can lie inside T grid but outside U, V, or f grid 
    8787   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   uo_e, vo_e 
    88    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ff_e, tt_e, fr_e 
     88   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   tt_e, fr_e 
    8989   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ua_e, va_e 
    9090   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ssh_e 
    9191   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   tmask_e, umask_e, vmask_e 
     92   REAl(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   rlon_e, rlat_e, ff_e 
     93      ! prepration to Nacho work 
     94      ! tt3d_e 
     95      ! uo3d_e 
     96      ! vo3d_e 
     97      ! 
    9298#if defined key_si3 || defined key_cice 
    9399   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   hi_e, ui_e, vi_e 
     
    181187         &      hi_e(0:jpi+1,0:jpj+1) ,                            & 
    182188#endif 
    183          &      ff_e(0:jpi+1,0:jpj+1) , fr_e(0:jpi+1,0:jpj+1)  ,   & 
     189         &      fr_e(0:jpi+1,0:jpj+1) ,                            & 
    184190         &      tt_e(0:jpi+1,0:jpj+1) , ssh_e(0:jpi+1,0:jpj+1) ,   & 
    185191         &      first_width(nclasses) , first_length(nclasses) ,   & 
     
    188194      icb_alloc = icb_alloc + ill 
    189195 
     196      ! prepration to Nacho work 
     197      ! tt3d_e 
     198      ! uo3d_e 
     199      ! vo3d_e 
     200      ! 
    190201      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), & 
    191          &      STAT=ill) 
     202         &      rlon_e(0:jpi+1,0:jpj+1) , rlat_e(0:jpi+1,0:jpj+1) , ff_e(0:jpi+1,0:jpj+1)   , STAT=ill) 
    192203      icb_alloc = icb_alloc + ill 
    193204 
  • NEMO/branches/2020/tickets_2494_2375/src/OCE/ICB/icbdyn.F90

    r10570 r13265  
    1414   USE dom_oce        ! NEMO ocean domain 
    1515   USE phycst         ! NEMO physical constants 
     16   USE in_out_manager                      ! IO parameters 
    1617   ! 
    1718   USE icb_oce        ! define iceberg arrays 
     
    155156         pt%xi   = zxi_n 
    156157         pt%yj   = zyj_n 
    157  
    158          ! update actual position 
    159          pt%lon  = icb_utl_bilin_x(glamt, pt%xi, pt%yj ) 
    160          pt%lat  = icb_utl_bilin(gphit, pt%xi, pt%yj, 'T' ) 
    161158 
    162159         berg => berg%next                         ! switch to the next berg 
     
    270267      ! Interpolate gridded fields to berg 
    271268      nknberg = berg%number(1) 
    272       CALL icb_utl_interp( pxi, pe1, zuo, zui, zua, zssh_x,                     & 
    273          &                 pyj, pe2, zvo, zvi, zva, zssh_y, zsst, zcn, zhi, zff ) 
     269      CALL icb_utl_interp( pxi, pyj, pe1=pe1, pe2=pe2,   &   ! scale factor 
     270         &                 puo=zuo, pui=zui, pua=zua,    &   ! oce/ice/atm velocities 
     271         &                 pvo=zvo, pvi=zvi, pva=zva,    &   ! oce/ice/atm velocities 
     272         &                 pssh_i=zssh_x, pssh_j=zssh_y, &   ! ssh gradient 
     273         &                 phi=zhi, pff=zff )               ! ice thickness and coriolis 
     274 
     275      IF (lwp) THEN 
     276         WRITE(numout,*) 'icbdyn ', pxi, pyj  
     277         WRITE(numout,*) pe1, zuo, zui, zua, zssh_x 
     278         WRITE(numout,*) pe2, zvo, zvi, zva, zssh_y 
     279         WRITE(numout,*) zhi, zff  
     280         WRITE(numout,*) '' 
     281      END IF 
     282 
    274283 
    275284      zM = berg%current_point%mass 
  • NEMO/branches/2020/tickets_2494_2375/src/OCE/ICB/icbini.F90

    r12489 r13265  
    8282      ff_e(:,:) = 0._wp   ;   tt_e(:,:) = 0._wp   ; 
    8383      fr_e(:,:) = 0._wp   ; 
     84      ! prepration to Nacho work 
     85      ! e3t_e 
     86      ! tt3d_e 
     87      ! uo3d_e 
     88      ! vo3d_e 
     89      ! 
    8490#if defined key_si3 
    8591      hi_e(:,:) = 0._wp   ; 
     
    242248      CALL lbc_lnk_icb( 'icbini', umask_e, 'T', +1._wp, 1, 1 ) 
    243249      CALL lbc_lnk_icb( 'icbini', vmask_e, 'T', +1._wp, 1, 1 ) 
    244       ! 
     250 
     251      ! definition of extended lat/lon array needed by icb_bilin_h 
     252      rlon_e(:,:) = 0._wp     ;  rlon_e(1:jpi,1:jpj) = glamt(:,:)  
     253      rlat_e(:,:) = 0._wp     ;  rlat_e(1:jpi,1:jpj) = gphit(:,:) 
     254      CALL lbc_lnk_icb( 'icbini', rlon_e, 'T', +1._wp, 1, 1 ) 
     255      CALL lbc_lnk_icb( 'icbini', rlat_e, 'T', +1._wp, 1, 1 ) 
     256      ! 
     257      ! definnitionn of extennded ff_f array needed by icb_utl_interp 
     258      ff_e(:,:) = 0._wp       ;  ff_e(1:jpi,1:jpj) = ff_f(:,:) 
     259      CALL lbc_lnk_icb( 'icbini', ff_e, 'F', +1._wp, 1, 1 ) 
     260 
    245261      ! assign each new iceberg with a unique number constructed from the processor number 
    246262      ! and incremented by the total number of processors 
     
    338354               localpt%xi = REAL( mig(ji), wp ) 
    339355               localpt%yj = REAL( mjg(jj), wp ) 
    340                localpt%lon = icb_utl_bilin(glamt, localpt%xi, localpt%yj, 'T' ) 
    341                localpt%lat = icb_utl_bilin(gphit, localpt%xi, localpt%yj, 'T' ) 
     356               CALL icb_utl_interp( localpt%xi, localpt%yj, plat=localpt%lat, plon=localpt%lon )    
    342357               localpt%mass      = rn_initial_mass     (iberg) 
    343358               localpt%thickness = rn_initial_thickness(iberg) 
  • NEMO/branches/2020/tickets_2494_2375/src/OCE/ICB/icbthm.F90

    r13226 r13265  
    8383         pt => this%current_point 
    8484         nknberg = this%number(1) 
    85          CALL icb_utl_interp( pt%xi, pt%e1, pt%uo, pt%ui, pt%ua, pt%ssh_x,   & 
    86             &                 pt%yj, pt%e2, pt%vo, pt%vi, pt%va, pt%ssh_y,   & 
    87             &                 pt%sst, pt%cn, pt%hi, zff ) 
     85         IF ( nn_sample_rate > 0 .AND. MOD(kt-1,nn_sample_rate) == 0 ) THEN 
     86            CALL icb_utl_interp( pt%xi, pt%yj, pe1=pt%e1, pe2=pt%e2,                 & 
     87               &                 puo=pt%uo, pui=pt%ui, pua=pt%ua, pssh_i=pt%ssh_x,   & 
     88               &                 pvo=pt%vo, pvi=pt%vi, pva=pt%va, pssh_j=pt%ssh_y,   & 
     89               &                 psst=pt%sst, pcn=pt%cn, phi=pt%hi,                  & 
     90               &                 plat=pt%lat, plon=pt%lon                            ) 
     91         ELSE 
     92            CALL icb_utl_interp( pt%xi, pt%yj,           &   ! position 
     93               &                 puo=pt%uo, pua=pt%ua,   &   ! oce/atm velocities 
     94               &                 pvo=pt%vo, pva=pt%va,   &   ! oce/atm velocities 
     95               &                 psst=pt%sst, pcn=pt%cn )    ! sst and ice concentration 
     96               ! preparation Nacho add t3d and uo, vo, to basal 
     97         END IF 
    8898         ! 
    8999         zSST = pt%sst 
  • NEMO/branches/2020/tickets_2494_2375/src/OCE/ICB/icbutl.F90

    r10702 r13265  
    88   !!            -    !                            Removal of mapping from another grid 
    99   !!            -    !  2011-04  (Alderson)       Split into separate modules 
     10   !!           4.2   !  2020-07  (P. Mathiot)     simplification of interpolation routine 
     11   !!                 !                            and add Nacho Merino work 
    1012   !!---------------------------------------------------------------------- 
    1113 
    1214   !!---------------------------------------------------------------------- 
    1315   !!   icb_utl_interp   : 
    14    !!   icb_utl_bilin    : 
    15    !!   icb_utl_bilin_e  : 
     16   !!   icb_utl_pos      : compute bottom left corner indice, weight and mask 
     17   !!   icb_utl_bilin_h  : interpolation field to icb position 
     18   !!   icb_utl_bilin_e  : interpolation of scale factor to icb position 
    1619   !!---------------------------------------------------------------------- 
    1720   USE par_oce                             ! ocean parameters 
     
    3336   PUBLIC   icb_utl_copy          ! routine called in icbstp module 
    3437   PUBLIC   icb_utl_interp        ! routine called in icbdyn, icbthm modules 
    35    PUBLIC   icb_utl_bilin         ! routine called in icbini, icbdyn modules 
    36    PUBLIC   icb_utl_bilin_x       ! routine called in icbdyn module 
     38   PUBLIC   icb_utl_bilin_h       ! routine called in icbdyn module 
    3739   PUBLIC   icb_utl_add           ! routine called in icbini.F90, icbclv, icblbc and icbrst modules 
    3840   PUBLIC   icb_utl_delete        ! routine called in icblbc, icbthm modules 
     
    7274      uo_e(1:jpi,1:jpj) = ssu_m(:,:) * umask(:,:,1) 
    7375      vo_e(1:jpi,1:jpj) = ssv_m(:,:) * vmask(:,:,1) 
    74       ff_e(1:jpi,1:jpj) = ff_f (:,:)  
    7576      tt_e(1:jpi,1:jpj) = sst_m(:,:) 
    7677      fr_e(1:jpi,1:jpj) = fr_i (:,:) 
     
    8081      CALL lbc_lnk_icb( 'icbutl', uo_e, 'U', -1._wp, 1, 1 ) 
    8182      CALL lbc_lnk_icb( 'icbutl', vo_e, 'V', -1._wp, 1, 1 ) 
    82       CALL lbc_lnk_icb( 'icbutl', ff_e, 'F', +1._wp, 1, 1 ) 
    8383      CALL lbc_lnk_icb( 'icbutl', ua_e, 'U', -1._wp, 1, 1 ) 
    8484      CALL lbc_lnk_icb( 'icbutl', va_e, 'V', -1._wp, 1, 1 ) 
    85       CALL lbc_lnk_icb( 'icbutl', fr_e, 'T', +1._wp, 1, 1 ) 
    86       CALL lbc_lnk_icb( 'icbutl', tt_e, 'T', +1._wp, 1, 1 ) 
    8785#if defined key_si3 
    8886      hi_e(1:jpi, 1:jpj) = hm_i (:,:)   
     
    9492      ssh_e(1:jpi, 1:jpj) = zssh_lead_m(:,:) * tmask(:,:,1) 
    9593      ! 
    96       CALL lbc_lnk_icb( 'icbutl', hi_e , 'T', +1._wp, 1, 1 ) 
    9794      CALL lbc_lnk_icb( 'icbutl', ui_e , 'U', -1._wp, 1, 1 ) 
    9895      CALL lbc_lnk_icb( 'icbutl', vi_e , 'V', -1._wp, 1, 1 ) 
    9996#else 
    100       ssh_e(1:jpi, 1:jpj) = ssh_m(:,:) * tmask(:,:,1) 
     97      ssh_e(1:jpi, 1:jpj) = ssh_m(:,:) * tmask(:,:,1)          
    10198#endif 
    102       CALL lbc_lnk_icb( 'icbutl', ssh_e, 'T', +1._wp, 1, 1 ) 
    10399      ! 
    104100   END SUBROUTINE icb_utl_copy 
    105101 
    106102 
    107    SUBROUTINE icb_utl_interp( pi, pe1, puo, pui, pua, pssh_i,   & 
    108       &                       pj, pe2, pvo, pvi, pva, pssh_j,   & 
    109       &                       psst, pcn, phi, pff            ) 
     103   SUBROUTINE icb_utl_interp( pi, pj, pe1, puo, pui, pua, pssh_i,   & 
     104      &                               pe2, pvo, pvi, pva, pssh_j,   & 
     105      &                               psst, pcn, phi, pff, plon, plat) 
    110106      !!---------------------------------------------------------------------- 
    111107      !!                  ***  ROUTINE icb_utl_interp  *** 
     
    125121      !!---------------------------------------------------------------------- 
    126122      REAL(wp), INTENT(in   ) ::   pi , pj                        ! position in (i,j) referential 
    127       REAL(wp), INTENT(  out) ::   pe1, pe2                       ! i- and j scale factors 
    128       REAL(wp), INTENT(  out) ::   puo, pvo, pui, pvi, pua, pva   ! ocean, ice and wind speeds 
    129       REAL(wp), INTENT(  out) ::   pssh_i, pssh_j                 ! ssh i- & j-gradients 
    130       REAL(wp), INTENT(  out) ::   psst, pcn, phi, pff            ! SST, ice concentration, ice thickness, Coriolis 
    131       ! 
     123      REAL(wp), INTENT(  out), OPTIONAL ::   pe1, pe2                       ! i- and j scale factors 
     124      REAL(wp), INTENT(  out), OPTIONAL ::   puo, pvo, pui, pvi, pua, pva   ! ocean, ice and wind speeds 
     125      REAL(wp), INTENT(  out), OPTIONAL ::   pssh_i, pssh_j                 ! ssh i- & j-gradients 
     126      REAL(wp), INTENT(  out), OPTIONAL ::   psst, pcn, phi, pff            ! SST, ice concentration, ice thickness, Coriolis 
     127      REAL(wp), INTENT(  out), OPTIONAL ::   plat, plon                     ! position 
     128      ! 
     129      REAL(wp), DIMENSION(4) :: zwT  , zwU  , zwV  , zwF   ! interpolation weight 
     130      REAL(wp), DIMENSION(4) :: zmskF, zmskU, zmskV, zmskT ! mask 
     131      REAL(wp), DIMENSION(4) :: zwTp, zmskTp, zwTm, zmskTm 
     132      INTEGER                :: iiT, iiU, iiV, iiF, ijT, ijU, ijV, ijF ! bottom left corner 
     133      INTEGER                :: iiTp, iiTm, ijTp, ijTm 
    132134      REAL(wp) ::   zcd, zmod       ! local scalars 
    133135      !!---------------------------------------------------------------------- 
    134  
    135       pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj )      ! scale factors 
    136       pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj ) 
    137       ! 
    138       puo  = icb_utl_bilin_h( uo_e, pi, pj, 'U', .false.  )    ! ocean velocities 
    139       pvo  = icb_utl_bilin_h( vo_e, pi, pj, 'V', .false.  ) 
    140       psst = icb_utl_bilin_h( tt_e, pi, pj, 'T', .true.   )    ! SST 
    141       pcn  = icb_utl_bilin_h( fr_e, pi, pj, 'T', .true.   )    ! ice concentration 
    142       pff  = icb_utl_bilin_h( ff_e, pi, pj, 'F', .false.  )    ! Coriolis parameter 
    143       ! 
    144       pua  = icb_utl_bilin_h( ua_e, pi, pj, 'U', .true.   )    ! 10m wind 
    145       pva  = icb_utl_bilin_h( va_e, pi, pj, 'V', .true.   )    ! here (ua,va) are stress => rough conversion from stress to speed 
    146       zcd  = 1.22_wp * 1.5e-3_wp                               ! air density * drag coefficient  
    147       zmod = 1._wp / MAX(  1.e-20, SQRT(  zcd * SQRT( pua*pua + pva*pva)  )  ) 
    148       pua  = pua * zmod                                       ! note: stress module=0 necessarly implies ua=va=0 
    149       pva  = pva * zmod 
    150  
     136      ! 
     137      ! get position, weight and mask  
     138      CALL icb_utl_pos( pi, pj, 'T', iiT, ijT, zwT, zmskT ) 
     139      CALL icb_utl_pos( pi, pj, 'U', iiU, ijU, zwU, zmskU ) 
     140      CALL icb_utl_pos( pi, pj, 'V', iiV, ijV, zwV, zmskV ) 
     141      CALL icb_utl_pos( pi, pj, 'F', iiF, ijF, zwF, zmskF ) 
     142      ! 
     143      ! metrics and coordinates 
     144      IF ( PRESENT(pe1 ) ) pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj )      ! scale factors 
     145      IF ( PRESENT(pe2 ) ) pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj ) 
     146      IF ( PRESENT(plon) ) plon= icb_utl_bilin_h( rlon_e, iiT, ijT, zwT, .true.  ) 
     147      IF ( PRESENT(plat) ) plat= icb_utl_bilin_h( rlat_e, iiT, ijT, zwT, .false. ) 
     148      ! 
     149      IF ( PRESENT(puo ) ) puo  = icb_utl_bilin_h( uo_e, iiU, ijU, zwU        , .false. )         ! ocean velocities 
     150      IF ( PRESENT(pvo ) ) pvo  = icb_utl_bilin_h( vo_e, iiV, ijV, zwV        , .false. )         ! 
     151      IF ( PRESENT(psst) ) psst = icb_utl_bilin_h( tt_e, iiT, ijT, zwT * zmskT, .false. ) ! sst 
     152      IF ( PRESENT(pcn ) ) pcn  = icb_utl_bilin_h( fr_e, iiT, ijT, zwT * zmskT, .false. ) ! ice concentration 
     153      IF ( PRESENT(pff ) ) pff  = icb_utl_bilin_h( ff_e, iiF, ijF, zwF        , .false. )         ! Coriolis parameter 
     154      ! 
     155      IF ( PRESENT(pua) .AND. PRESENT(pva) ) THEN 
     156         pua  = icb_utl_bilin_h( ua_e, iiU, ijU, zwU * zmskU, .false. ) ! 10m wind 
     157         pva  = icb_utl_bilin_h( va_e, iiV, ijV, zwV * zmskV, .false. ) ! here (ua,va) are stress => rough conversion from stress to speed 
     158         zcd  = 1.22_wp * 1.5e-3_wp                               ! air density * drag coefficient  
     159         zmod = 1._wp / MAX(  1.e-20, SQRT(  zcd * SQRT( pua*pua + pva*pva)  )  ) 
     160         pua  = pua * zmod                                       ! note: stress module=0 necessarly implies ua=va=0 
     161         pva  = pva * zmod 
     162      END IF 
     163      ! 
    151164#if defined key_si3 
    152       pui = icb_utl_bilin_h( ui_e , pi, pj, 'U', .false. )    ! sea-ice velocities 
    153       pvi = icb_utl_bilin_h( vi_e , pi, pj, 'V', .false. ) 
    154       phi = icb_utl_bilin_h( hi_e , pi, pj, 'T', .true.  )    ! ice thickness 
     165      IF ( PRESENT(pui) ) pui = icb_utl_bilin_h( ui_e , iiU, ijU, zwU        , .false. ) ! sea-ice velocities 
     166      IF ( PRESENT(pvi) ) pvi = icb_utl_bilin_h( vi_e , iiV, ijV, zwV        , .false. ) 
     167      IF ( PRESENT(phi) ) phi = icb_utl_bilin_h( hi_e , iiT, ijT, zwT * zmskT, .false. ) ! ice thickness 
    155168#else 
    156       pui = 0._wp 
    157       pvi = 0._wp 
    158       phi = 0._wp 
     169      IF ( PRESENT(pui) ) pui = 0._wp 
     170      IF ( PRESENT(pvi) ) pvi = 0._wp 
     171      IF ( PRESENT(phi) ) phi = 0._wp 
    159172#endif 
    160  
     173      ! 
    161174      ! Estimate SSH gradient in i- and j-direction (centred evaluation) 
    162       pssh_i = ( icb_utl_bilin_h( ssh_e, pi+0.1_wp, pj, 'T', .true. ) -   & 
    163          &       icb_utl_bilin_h( ssh_e, pi-0.1_wp, pj, 'T', .true. )  ) / ( 0.2_wp * pe1 ) 
    164       pssh_j = ( icb_utl_bilin_h( ssh_e, pi, pj+0.1_wp, 'T', .true. ) -   & 
    165          &       icb_utl_bilin_h( ssh_e, pi, pj-0.1_wp, 'T', .true. )  ) / ( 0.2_wp * pe2 ) 
     175      IF ( PRESENT(pssh_i) .AND. PRESENT(pssh_j) ) THEN 
     176         CALL icb_utl_pos( pi+0.1, pj    , 'T', iiTp, ijTp, zwTp, zmskTp ) 
     177         CALL icb_utl_pos( pi-0.1, pj    , 'T', iiTm, ijTm, zwTm, zmskTm ) 
     178         ! 
     179         IF ( .NOT. PRESENT(pe1) ) pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj ) 
     180         pssh_i = ( icb_utl_bilin_h( ssh_e, iiTp, ijTp, zwTp*zmskTp, .false. ) -   & 
     181            &       icb_utl_bilin_h( ssh_e, iiTm, ijTm, zwTm*zmskTm, .false. )  ) / ( 0.2_wp * pe1 ) 
     182         ! 
     183         CALL icb_utl_pos( pi    , pj+0.1, 'T', iiTp, ijTp, zwTp, zmskTp ) 
     184         CALL icb_utl_pos( pi    , pj-0.1, 'T', iiTm, ijTm, zwTm, zmskTm ) 
     185         ! 
     186         IF ( .NOT. PRESENT(pe2) ) pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj ) 
     187         pssh_j = ( icb_utl_bilin_h( ssh_e, iiTp, ijTp, zwTp*zmskTp, .false. ) -   & 
     188            &       icb_utl_bilin_h( ssh_e, iiTm, ijTm, zwTm*zmskTm, .false. )  ) / ( 0.2_wp * pe2 ) 
     189      END IF 
    166190      ! 
    167191   END SUBROUTINE icb_utl_interp 
    168192 
    169  
    170    REAL(wp) FUNCTION icb_utl_bilin_h( pfld, pi, pj, cd_type, plmask ) 
     193   SUBROUTINE icb_utl_pos( pi, pj, cd_type, kii, kij, pw, pmsk ) 
    171194      !!---------------------------------------------------------------------- 
    172195      !!                  ***  FUNCTION icb_utl_bilin  *** 
     
    179202      !! 
    180203      !!---------------------------------------------------------------------- 
    181       REAL(wp), DIMENSION(0:jpi+1,0:jpj+1), INTENT(in) ::   pfld      ! field to be interpolated 
    182       REAL(wp)                            , INTENT(in) ::   pi, pj    ! targeted coordinates in (i,j) referential 
    183       CHARACTER(len=1)                    , INTENT(in) ::   cd_type   ! type of pfld array grid-points: = T , U , V or F points 
    184       LOGICAL                             , INTENT(in) ::   plmask    ! special treatment of mask point 
    185       ! 
    186       INTEGER  ::   ii, ij   ! local integer 
    187       REAL(wp) ::   zi, zj   ! local real 
    188       REAL(wp) :: zw1, zw2, zw3, zw4 
    189       REAL(wp), DIMENSION(4) :: zmask 
     204      REAL(wp)              , INTENT(IN)  ::   pi, pj    ! targeted coordinates in (i,j) referential 
     205      CHARACTER(len=1)      , INTENT(IN)  ::   cd_type   ! point type 
     206      REAL(wp), DIMENSION(4), INTENT(OUT) ::   pw, pmsk  ! weight and mask 
     207      INTEGER ,               INTENT(OUT) ::   kii, kij  ! bottom left corner position in local domain 
     208      ! 
     209      REAL(wp) :: zwi, zwj ! distance to bottom left corner 
     210      INTEGER  :: ierr  
     211      ! 
    190212      !!---------------------------------------------------------------------- 
    191213      ! 
     
    195217         ! since we're looking for four T points containing quadrant we're in of  
    196218         ! current T cell 
    197          ii = MAX(0, INT( pi     )) 
    198          ij = MAX(0, INT( pj     ))    ! T-point 
    199          zi = pi - REAL(ii,wp) 
    200          zj = pj - REAL(ij,wp) 
     219         kii = MAX(0, INT( pi     )) 
     220         kij = MAX(0, INT( pj     ))    ! T-point 
     221         zwi = pi - REAL(kii,wp) 
     222         zwj = pj - REAL(kij,wp) 
    201223      CASE ( 'U' ) 
    202          ii = MAX(0, INT( pi-0.5_wp )) 
    203          ij = MAX(0, INT( pj     ))    ! U-point 
    204          zi = pi - 0.5_wp - REAL(ii,wp) 
    205          zj = pj - REAL(ij,wp) 
     224         kii = MAX(0, INT( pi-0.5_wp )) 
     225         kij = MAX(0, INT( pj     ))    ! U-point 
     226         zwi = pi - 0.5_wp - REAL(kii,wp) 
     227         zwj = pj - REAL(kij,wp) 
    206228      CASE ( 'V' ) 
    207          ii = MAX(0, INT( pi     )) 
    208          ij = MAX(0, INT( pj-0.5_wp ))    ! V-point 
    209          zi = pi - REAL(ii,wp) 
    210          zj = pj - 0.5_wp - REAL(ij,wp) 
     229         kii = MAX(0, INT( pi     )) 
     230         kij = MAX(0, INT( pj-0.5_wp ))    ! V-point 
     231         zwi = pi - REAL(kii,wp) 
     232         zwj = pj - 0.5_wp - REAL(kij,wp) 
    211233      CASE ( 'F' ) 
    212          ii = MAX(0, INT( pi-0.5_wp )) 
    213          ij = MAX(0, INT( pj-0.5_wp ))    ! F-point 
    214          zi = pi - 0.5_wp - REAL(ii,wp) 
    215          zj = pj - 0.5_wp - REAL(ij,wp) 
     234         kii = MAX(0, INT( pi-0.5_wp )) 
     235         kij = MAX(0, INT( pj-0.5_wp ))    ! F-point 
     236         zwi = pi - 0.5_wp - REAL(kii,wp) 
     237         zwj = pj - 0.5_wp - REAL(kij,wp) 
    216238      END SELECT 
     239      ! 
     240      ! compute weight 
     241      pw(1) = (1._wp-zwi) * (1._wp-zwj) 
     242      pw(2) =        zwi  * (1._wp-zwj) 
     243      pw(3) = (1._wp-zwi) *        zwj 
     244      pw(4) =        zwi  *        zwj 
     245      ! 
     246      ! find position in this processor. Prevent near edge problems (see #1389) 
     247      ! 
     248      IF (TRIM(cd_type) == 'T' ) THEN 
     249         ierr = 0 
     250         IF    ( kii <  mig( 1 ) ) THEN   ;  ierr = ierr + 1 
     251         ELSEIF( kii >= mig(jpi) ) THEN   ;  ierr = ierr + 1 
     252         ENDIF 
     253         ! 
     254         IF    ( kij <  mjg( 1 ) ) THEN   ;   ierr = ierr + 1 
     255         ELSEIF( kij >= mjg(jpj) ) THEN   ;   ierr = ierr + 1 
     256         ENDIF 
     257         ! 
     258         IF ( ierr > 0 ) THEN 
     259            WRITE(numout,*) 'bottom left corner T point out of bound' 
     260            WRITE(numout,*) kii, mig( 1 ), mig(jpi) 
     261            WRITE(numout,*) kij, mjg( 1 ), mjg(jpj) 
     262            CALL ctl_stop('STOP','icb_utl_bilin_h: an icebergs coordinates is out of valid range (out of bound error)') 
     263         END IF 
     264      END IF 
    217265      ! 
    218266      ! find position in this processor. Prevent near edge problems (see #1389) 
    219267      ! (PM) will be useless if extra halo is used in NEMO 
    220268      ! 
    221       IF    ( ii <= mig(1)-1 ) THEN   ;   ii = 0 
    222       ELSEIF( ii  > mig(jpi) ) THEN   ;   ii = jpi 
    223       ELSE                            ;   ii = mi1(ii) 
    224       ENDIF 
    225       IF    ( ij <= mjg(1)-1 ) THEN   ;   ij = 0 
    226       ELSEIF( ij  > mjg(jpj) ) THEN   ;   ij = jpj 
    227       ELSE                            ;   ij = mj1(ij) 
     269      IF    ( kii <= mig(1)-1 ) THEN   ;   kii = 0 
     270      ELSEIF( kii  > mig(jpi) ) THEN   ;   kii = jpi 
     271      ELSE                             ;   kii = mi1(kii) 
     272      ENDIF 
     273      IF    ( kij <= mjg(1)-1 ) THEN   ;   kij = 0 
     274      ELSEIF( kij  > mjg(jpj) ) THEN   ;   kij = jpj 
     275      ELSE                             ;   kij = mj1(kij) 
    228276      ENDIF 
    229277      ! 
    230278      ! define mask array  
    231       IF (plmask) THEN 
    232          ! land value is not used in the interpolation 
    233          SELECT CASE ( cd_type ) 
    234          CASE ( 'T' ) 
    235             zmask = (/tmask_e(ii,ij), tmask_e(ii+1,ij), tmask_e(ii,ij+1), tmask_e(ii+1,ij+1)/) 
    236          CASE ( 'U' ) 
    237             zmask = (/umask_e(ii,ij), umask_e(ii+1,ij), umask_e(ii,ij+1), umask_e(ii+1,ij+1)/) 
    238          CASE ( 'V' ) 
    239             zmask = (/vmask_e(ii,ij), vmask_e(ii+1,ij), vmask_e(ii,ij+1), vmask_e(ii+1,ij+1)/) 
    240          CASE ( 'F' ) 
    241             ! F case only used for coriolis, ff_f is not mask so zmask = 1 
    242             zmask = 1. 
    243          END SELECT 
    244       ELSE 
    245          ! land value is used during interpolation 
    246          zmask = 1. 
    247       END iF 
    248       ! 
    249       ! compute weight 
    250       zw1 = zmask(1) * (1._wp-zi) * (1._wp-zj) 
    251       zw2 = zmask(2) *        zi  * (1._wp-zj) 
    252       zw3 = zmask(3) * (1._wp-zi) *        zj 
    253       zw4 = zmask(4) *        zi  *        zj 
    254       ! 
    255       ! compute interpolated value 
    256       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)  
    257       ! 
    258    END FUNCTION icb_utl_bilin_h 
    259  
    260  
    261    REAL(wp) FUNCTION icb_utl_bilin( pfld, pi, pj, cd_type ) 
     279      ! land value is not used in the interpolation 
     280      SELECT CASE ( cd_type ) 
     281      CASE ( 'T' ) 
     282         pmsk = (/tmask_e(kii,kij), tmask_e(kii+1,kij), tmask_e(kii,kij+1), tmask_e(kii+1,kij+1)/) 
     283      CASE ( 'U' ) 
     284         pmsk = (/umask_e(kii,kij), umask_e(kii+1,kij), umask_e(kii,kij+1), umask_e(kii+1,kij+1)/) 
     285      CASE ( 'V' ) 
     286         pmsk = (/vmask_e(kii,kij), vmask_e(kii+1,kij), vmask_e(kii,kij+1), vmask_e(kii+1,kij+1)/) 
     287      CASE ( 'F' ) 
     288         ! F case only used for coriolis, ff_f is not mask so zmask = 1 
     289         pmsk = 1. 
     290      END SELECT 
     291   END SUBROUTINE icb_utl_pos 
     292 
     293   REAL(wp) FUNCTION icb_utl_bilin_h( pfld, pii, pij, pw, pllon ) 
    262294      !!---------------------------------------------------------------------- 
    263295      !!                  ***  FUNCTION icb_utl_bilin  *** 
    264296      !! 
    265297      !! ** Purpose :   bilinear interpolation at berg location depending on the grid-point type 
     298      !!                this version deals with extra halo points 
    266299      !! 
    267300      !!       !!gm  CAUTION an optional argument should be added to handle 
     
    269302      !! 
    270303      !!---------------------------------------------------------------------- 
    271       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pfld      ! field to be interpolated 
    272       REAL(wp)                    , INTENT(in) ::   pi, pj    ! targeted coordinates in (i,j) referential 
    273       CHARACTER(len=1)            , INTENT(in) ::   cd_type   ! type of pfld array grid-points: = T , U , V or F points 
    274       ! 
    275       INTEGER  ::   ii, ij   ! local integer 
    276       REAL(wp) ::   zi, zj   ! local real 
    277       !!---------------------------------------------------------------------- 
    278       ! 
    279       SELECT CASE ( cd_type ) 
    280          CASE ( 'T' ) 
    281             ! note that here there is no +0.5 added 
    282             ! since we're looking for four T points containing quadrant we're in of  
    283             ! current T cell 
    284             ii = MAX(1, INT( pi     )) 
    285             ij = MAX(1, INT( pj     ))    ! T-point 
    286             zi = pi - REAL(ii,wp) 
    287             zj = pj - REAL(ij,wp) 
    288          CASE ( 'U' ) 
    289             ii = MAX(1, INT( pi-0.5 )) 
    290             ij = MAX(1, INT( pj     ))    ! U-point 
    291             zi = pi - 0.5 - REAL(ii,wp) 
    292             zj = pj - REAL(ij,wp) 
    293          CASE ( 'V' ) 
    294             ii = MAX(1, INT( pi     )) 
    295             ij = MAX(1, INT( pj-0.5 ))    ! V-point 
    296             zi = pi - REAL(ii,wp) 
    297             zj = pj - 0.5 - REAL(ij,wp) 
    298          CASE ( 'F' ) 
    299             ii = MAX(1, INT( pi-0.5 )) 
    300             ij = MAX(1, INT( pj-0.5 ))    ! F-point 
    301             zi = pi - 0.5 - REAL(ii,wp) 
    302             zj = pj - 0.5 - REAL(ij,wp) 
    303       END SELECT 
    304       ! 
    305       ! find position in this processor. Prevent near edge problems (see #1389) 
    306       IF    ( ii < mig( 1 ) ) THEN   ;   ii = 1 
    307       ELSEIF( ii > mig(jpi) ) THEN   ;   ii = jpi 
    308       ELSE                           ;   ii = mi1(ii) 
    309       ENDIF 
    310       IF    ( ij < mjg( 1 ) ) THEN   ;   ij = 1 
    311       ELSEIF( ij > mjg(jpj) ) THEN   ;   ij = jpj 
    312       ELSE                           ;   ij  = mj1(ij) 
    313       ENDIF 
    314       ! 
    315       IF( ii == jpi )   ii = ii-1       
    316       IF( ij == jpj )   ij = ij-1 
    317       ! 
    318       icb_utl_bilin = ( pfld(ii,ij  ) * (1.-zi) + pfld(ii+1,ij  ) * zi ) * (1.-zj)   & 
    319          &          + ( pfld(ii,ij+1) * (1.-zi) + pfld(ii+1,ij+1) * zi ) *     zj 
    320       ! 
    321    END FUNCTION icb_utl_bilin 
    322  
    323  
    324    REAL(wp) FUNCTION icb_utl_bilin_x( pfld, pi, pj ) 
    325       !!---------------------------------------------------------------------- 
    326       !!                  ***  FUNCTION icb_utl_bilin_x  *** 
    327       !! 
    328       !! ** Purpose :   bilinear interpolation at berg location depending on the grid-point type 
    329       !!                Special case for interpolating longitude 
    330       !! 
    331       !!       !!gm  CAUTION an optional argument should be added to handle 
    332       !!             the slip/no-slip conditions  ==>>> to be done later 
    333       !! 
    334       !!---------------------------------------------------------------------- 
    335       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pfld      ! field to be interpolated 
    336       REAL(wp)                    , INTENT(in) ::   pi, pj    ! targeted coordinates in (i,j) referential 
    337       ! 
    338       INTEGER                                  ::   ii, ij   ! local integer 
    339       REAL(wp)                                 ::   zi, zj   ! local real 
    340       REAL(wp)                                 ::   zret     ! local real 
    341       REAL(wp), DIMENSION(4)                   ::   z4 
    342       !!---------------------------------------------------------------------- 
    343       ! 
    344       ! note that here there is no +0.5 added 
    345       ! since we're looking for four T points containing quadrant we're in of  
    346       ! current T cell 
    347       ii = MAX(1, INT( pi     )) 
    348       ij = MAX(1, INT( pj     ))    ! T-point 
    349       zi = pi - REAL(ii,wp) 
    350       zj = pj - REAL(ij,wp) 
    351       ! 
    352       ! find position in this processor. Prevent near edge problems (see #1389) 
    353       IF    ( ii < mig( 1 ) ) THEN   ;   ii = 1 
    354       ELSEIF( ii > mig(jpi) ) THEN   ;   ii = jpi 
    355       ELSE                           ;   ii = mi1(ii) 
    356       ENDIF 
    357       IF    ( ij < mjg( 1 ) ) THEN   ;   ij = 1 
    358       ELSEIF( ij > mjg(jpj) ) THEN   ;   ij = jpj 
    359       ELSE                           ;   ij  = mj1(ij) 
    360       ENDIF 
    361       ! 
    362       IF( ii == jpi )   ii = ii-1       
    363       IF( ij == jpj )   ij = ij-1 
    364       ! 
    365       z4(1) = pfld(ii  ,ij  ) 
    366       z4(2) = pfld(ii+1,ij  ) 
    367       z4(3) = pfld(ii  ,ij+1) 
    368       z4(4) = pfld(ii+1,ij+1) 
    369       IF( MAXVAL(z4) - MINVAL(z4) > 90._wp ) THEN 
    370          WHERE( z4 < 0._wp ) z4 = z4 + 360._wp 
    371       ENDIF 
    372       ! 
    373       zret = (z4(1) * (1.-zi) + z4(2) * zi) * (1.-zj) + (z4(3) * (1.-zi) + z4(4) * zi) * zj 
    374       IF( zret > 180._wp ) zret = zret - 360._wp 
    375       icb_utl_bilin_x = zret 
    376       ! 
    377    END FUNCTION icb_utl_bilin_x 
    378  
     304      REAL(wp), DIMENSION(0:jpi+1,0:jpj+1), INTENT(in) ::   pfld      ! field to be interpolated 
     305      REAL(wp), DIMENSION(4)              , INTENT(in) ::   pw        ! weight 
     306      LOGICAL                             , INTENT(in) ::   pllon     ! input data is a longitude 
     307      INTEGER ,                             INTENT(in) ::   pii, pij  ! bottom left corner 
     308      ! 
     309      REAL(wp), DIMENSION(4) :: zdat ! input data 
     310      !!---------------------------------------------------------------------- 
     311      ! 
     312      ! data 
     313      zdat(1) = pfld(pii  ,pij  ) 
     314      zdat(2) = pfld(pii+1,pij  ) 
     315      zdat(3) = pfld(pii  ,pij+1) 
     316      zdat(4) = pfld(pii+1,pij+1) 
     317      ! 
     318      IF( pllon .AND. MAXVAL(zdat) - MINVAL(zdat) > 90._wp ) THEN 
     319         WHERE( zdat < 0._wp ) zdat = zdat + 360._wp 
     320      ENDIF 
     321      ! 
     322      ! compute interpolated value 
     323      icb_utl_bilin_h = ( zdat(1)*pw(1) + zdat(2)*pw(2) + zdat(3)*pw(3) + zdat(4)*pw(4) ) / MAX(1.e-20, pw(1)+pw(2)+pw(3)+pw(4))  
     324      ! 
     325      IF( pllon .AND. icb_utl_bilin_h > 180._wp ) icb_utl_bilin_h = icb_utl_bilin_h - 360._wp 
     326      ! 
     327   END FUNCTION icb_utl_bilin_h 
    379328 
    380329   REAL(wp) FUNCTION icb_utl_bilin_e( pet, peu, pev, pef, pi, pj ) 
     
    387336      !!---------------------------------------------------------------------- 
    388337      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pet, peu, pev, pef   ! horizontal scale factor to be interpolated at t-,u-,v- & f-pts 
    389       REAL(wp)                , INTENT(in) ::   pi, pj               ! targeted coordinates in (i,j) referential 
    390       ! 
    391       INTEGER  ::   ii, ij, icase, ierr   ! local integer 
     338      REAL(wp)                , INTENT(IN) ::   pi , pj              ! iceberg position 
    392339      ! 
    393340      ! weights corresponding to corner points of a T cell quadrant 
    394341      REAL(wp) ::   zi, zj          ! local real 
     342      INTEGER  ::   ii, ij          ! bottom left corner coordinate in local domain 
    395343      ! 
    396344      ! values at corner points of a T cell quadrant 
     
    399347      !!---------------------------------------------------------------------- 
    400348      ! 
     349      ! cannot used iiT because need ii/ij reltaive to global indices not local one 
    401350      ii = MAX(1, INT( pi ))   ;   ij = MAX(1, INT( pj ))            ! left bottom T-point (i,j) indices 
    402  
     351      !  
    403352      ! fractional box spacing 
    404353      ! 0   <= zi < 0.5, 0   <= zj < 0.5   -->  NW quadrant of current T cell 
     
    410359      zj = pj - REAL(ij,wp) 
    411360 
    412       ! find position in this processor. Prevent near edge problems (see #1389) 
    413       ! 
    414       ierr = 0 
    415       IF    ( ii < mig( 1 ) ) THEN   ;   ii = 1       ; ierr = ierr + 1 
    416       ELSEIF( ii > mig(jpi) ) THEN   ;   ii = jpi     ; ierr = ierr + 1 
    417       ELSE                           ;   ii = mi1(ii) 
    418       ENDIF 
    419       IF    ( ij < mjg( 1 ) ) THEN   ;   ij = 1       ; ierr = ierr + 1 
    420       ELSEIF( ij > mjg(jpj) ) THEN   ;   ij = jpj     ; ierr = ierr + 1 
    421       ELSE                           ;   ij  = mj1(ij) 
    422       ENDIF 
    423       ! 
    424       IF( ii == jpi ) THEN ; ii = ii-1 ; ierr = ierr + 1 ; END IF      
    425       IF( ij == jpj ) THEN ; ij = ij-1 ; ierr = ierr + 1 ; END IF 
    426       ! 
    427       IF ( ierr > 0 ) CALL ctl_stop('STOP','icb_utl_bilin_e: an icebergs coordinates is out of valid range (out of bound error)') 
     361      ! conversion to local domain (no need to do a sanity check already done in icbpos) 
     362      ii = mi1(ii) 
     363      ij = mj1(ij) 
     364 
     365      IF (lwp) THEN 
     366         WRITE(numout,*) 'bilin_e ', pi, pj 
     367         WRITE(numout,*) ii, ij, zi, zj 
     368         WRITE(numout,*) '' 
     369      END IF 
     370 
    428371      ! 
    429372      IF(    0.0_wp <= zi .AND. zi < 0.5_wp   ) THEN 
Note: See TracChangeset for help on using the changeset viewer.