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 13374 for NEMO/branches – NEMO

Changeset 13374 for NEMO/branches


Ignore:
Timestamp:
2020-08-03T15:48:40+02:00 (4 years ago)
Author:
mathiot
Message:

ticket #1900: fix issue about mask management in icb_utl_bilin_3d_h; add option to ground icb if icb bottom lvl hit oce bottom lvl

Location:
NEMO/branches/2020/tickets_icb_1900
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/tickets_icb_1900/cfgs/SHARED/namelist_ref

    r13357 r13374  
    616616   ! 
    617617   ln_M2016                = .false. ! use Merino et al. (2016) modification (use of 3d ocean data instead of only sea surface data) 
     618      ln_icb_grd           = .false. ! ground icb when icb bottom level hit oce bottom level (need ln_M2016 to be activated) 
    618619   ! 
    619620   cn_dir      = './'      !  root directory for the calving data location 
  • NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icb_oce.F90

    r13357 r13374  
    129129   LOGICAL , PUBLIC ::   ln_time_average_weight          !: Time average the weight on the ocean    !!gm I don't understand that ! 
    130130   REAL(wp), PUBLIC ::   rn_speed_limit                  !: CFL speed limit for a berg 
    131    LOGICAL , PUBLIC ::   ln_M2016                        !: use Nacho's Merino 2016 work 
     131   LOGICAL , PUBLIC ::   ln_M2016, ln_icb_grd            !: use Nacho's Merino 2016 work 
    132132   ! 
    133133   ! restart 
  • NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbdyn.F90

    r13359 r13374  
    9898         zyj2 = zyj1 + zdt_2 * zv1          ;   zvvel2 = zvvel1 + zdt_2 * zay1 
    9999         ! 
    100          CALL icb_ground( zxi2, zxi1, zu1,   & 
    101             &             zyj2, zyj1, zv1, ll_bounced ) 
     100         CALL icb_ground( berg, zxi2, zxi1, zu1,   & 
     101            &                   zyj2, zyj1, zv1, ll_bounced ) 
    102102 
    103103         !                                         !**   A2 = A(X2,V2) 
     
    114114         zyj3  = zyj1  + zdt_2 * zv2   ;   zvvel3 = zvvel1 + zdt_2 * zay2 
    115115         ! 
    116          CALL icb_ground( zxi3, zxi1, zu3,   & 
    117             &                zyj3, zyj1, zv3, ll_bounced ) 
     116         CALL icb_ground( berg, zxi3, zxi1, zu3,   & 
     117            &                   zyj3, zyj1, zv3, ll_bounced ) 
    118118 
    119119         !                                         !**   A3 = A(X3,V3) 
     
    130130         zyj4 = zyj1 + zdt * zv3   ;   zvvel4 = zvvel1 + zdt * zay3 
    131131 
    132          CALL icb_ground( zxi4, zxi1, zu4,   & 
    133             &             zyj4, zyj1, zv4, ll_bounced ) 
     132         CALL icb_ground( berg, zxi4, zxi1, zu4,   & 
     133            &                   zyj4, zyj1, zv4, ll_bounced ) 
    134134 
    135135         !                                         !**   A4 = A(X4,V4) 
     
    149149         zvvel_n = pt%vvel + zdt_6 * (  zay1 + 2.*(zay2 + zay3) + zay4 ) 
    150150 
    151          CALL icb_ground( zxi_n, zxi1, zuvel_n,   & 
    152             &             zyj_n, zyj1, zvvel_n, ll_bounced ) 
     151         CALL icb_ground( berg, zxi_n, zxi1, zuvel_n,   & 
     152            &                   zyj_n, zyj1, zvvel_n, ll_bounced ) 
    153153 
    154154         pt%uvel = zuvel_n                        !** save in berg structure 
     
    164164 
    165165 
    166    SUBROUTINE icb_ground( pi, pi0, pu,   & 
    167       &                   pj, pj0, pv, ld_bounced ) 
     166   SUBROUTINE icb_ground( berg, pi, pi0, pu,   & 
     167      &                         pj, pj0, pv, ld_bounced ) 
    168168      !!---------------------------------------------------------------------- 
    169169      !!                  ***  ROUTINE icb_ground  *** 
     
    174174      !!                NB two possibilities available one of which is hard-coded here 
    175175      !!---------------------------------------------------------------------- 
     176      TYPE(iceberg ), POINTER, INTENT(in   ) ::   berg             ! berg 
     177      ! 
    176178      REAL(wp), INTENT(inout) ::   pi , pj      ! current iceberg position 
    177179      REAL(wp), INTENT(in   ) ::   pi0, pj0     ! previous iceberg position 
     
    181183      INTEGER  ::   ii, ii0 
    182184      INTEGER  ::   ij, ij0 
     185      INTEGER  ::   ikb 
    183186      INTEGER  ::   ibounce_method 
     187      ! 
     188      REAL(wp) :: zD  
     189      REAL(wp), DIMENSION(jpk) :: ze3t 
    184190      !!---------------------------------------------------------------------- 
    185191      ! 
     
    198204      ! 
    199205      ! assume icb is grounded if tmask(ii,ij,1) or tmask(ii,ij,ikb), depending of the option is not 0 
    200       !IF ( ln_icb_ground ) THEN 
    201       !   ! interpol needed data 
    202       !   CALL icb_utl_interp( pxi, pyj, pe3t=ze3t )   ! 3d velocities 
    203       !   
    204       !   !compute bottom level 
    205       !   CALL icb_utl_getkb( ikb, ze3t, zD ) 
    206       ! 
    207       !   IF(  tmask(ii,ij,ikb)  /=   0._wp  )   RETURN           ! berg reach a new t-cell, but an ocean one 
    208       !ELSE 
    209       IF(  tmask(ii,ij,1)  /=   0._wp  )   RETURN           ! berg reach a new t-cell, but an ocean one 
    210       !END IF 
     206      IF ( ln_M2016 .AND. ln_icb_grd ) THEN 
     207         ! 
     208         ! draught (keel depth) 
     209         zD = rho_berg_1_oce * berg%current_point%thickness 
     210         ! 
     211         ! interpol needed data 
     212         CALL icb_utl_interp( pi, pj, pe3t=ze3t ) 
     213         !  
     214         !compute bottom level 
     215         CALL icb_utl_getkb( ikb, ze3t, zD ) 
     216         ! 
     217         ! berg reach a new t-cell, but an ocean one 
     218         IF(  tmask(ii,ij,ikb) /= 0._wp .AND. tmask(ii,ij,1) /= 0._wp ) RETURN 
     219         ! 
     220      ELSE 
     221         IF(  tmask(ii,ij,1)  /=   0._wp  )   RETURN           ! berg reach a new t-cell, but an ocean one 
     222      END IF 
    211223      ! 
    212224      ! From here, berg have reach land: treat grounding/bouncing 
  • NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbini.F90

    r13365 r13374  
    415415         &              ln_time_average_weight          , nn_test_icebergs    , rn_test_box          ,   & 
    416416         &              ln_use_calving , rn_speed_limit , cn_dir, sn_icb      , ln_M2016             ,   & 
    417          &              cn_icbrst_indir, cn_icbrst_in   , cn_icbrst_outdir    , cn_icbrst_out 
     417         &              cn_icbrst_indir, cn_icbrst_in   , cn_icbrst_outdir    , cn_icbrst_out        ,   & 
     418         &              ln_icb_grd 
    418419      !!---------------------------------------------------------------------- 
    419420 
     
    494495            &                    'bits_erosion_fraction = ', rn_bits_erosion_fraction 
    495496 
     497         WRITE(numout,*) '   Use icb module modification from Merino et al. (2016) : ln_M2016 = ', ln_M2016 
     498         WRITE(numout,*) '       ground icebergs if icb bottom lvl hit the oce bottom level : ln_icb_grd = ', ln_icb_grd 
     499 
    496500         WRITE(numout,*) '   Shift of sea-ice concentration in erosion flux modulation ',   & 
    497501            &                    '(0<sicn_shift<1)    rn_sicn_shift  = ', rn_sicn_shift 
  • NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbutl.F90

    r13365 r13374  
    164164      REAL(wp), DIMENSION(4) :: zmskF, zmskU, zmskV, zmskT ! mask 
    165165      REAL(wp), DIMENSION(4) :: zwTp, zmskTp, zwTm, zmskTm 
     166      REAL(wp), DIMENSION(4,jpk) :: zw1d 
    166167      INTEGER                :: iiT, iiU, iiV, iiF, ijT, ijU, ijV, ijF ! bottom left corner 
    167168      INTEGER                :: iiTp, iiTm, ijTp, ijTm 
     
    225226      ! 3d interpolation 
    226227      IF ( PRESENT(puoce) .AND. PRESENT(pvoce) ) THEN 
    227          puoce(:) = icb_utl_bilin_h( uoce_e , iiU, ijU, zwU ) 
    228          pvoce(:) = icb_utl_bilin_h( voce_e , iiV, ijV, zwV ) 
     228         ! no need to mask as 0 is a valid data for land 
     229         zw1d(1,:) = zwU(1) ; zw1d(2,:) = zwU(2) ; zw1d(3,:) = zwU(3) ; zw1d(4,:) = zwU(4) ; 
     230         puoce(:) = icb_utl_bilin_h( uoce_e , iiU, ijU, zw1d ) 
     231 
     232         zw1d(1,:) = zwV(1) ; zw1d(2,:) = zwV(2) ; zw1d(3,:) = zwV(3) ; zw1d(4,:) = zwV(4) ; 
     233         pvoce(:) = icb_utl_bilin_h( voce_e , iiV, ijV, zw1d ) 
    229234      END IF 
    230       IF ( PRESENT(ptoce) ) ptoce(:) = icb_utl_bilin_h( toce_e , iiT, ijT, zwT * zmskT ) 
     235 
     236      IF ( PRESENT(ptoce) ) THEN 
     237         ! for temperature we need to mask the weight properly 
     238         ! no need of extra halo as it is a T point variable 
     239         zw1d(1,:) = tmask(iiT  ,ijT  ,:) * zwT(1) * zmskT(1) 
     240         zw1d(2,:) = tmask(iiT+1,ijT  ,:) * zwT(2) * zmskT(2) 
     241         zw1d(3,:) = tmask(iiT  ,ijT+1,:) * zwT(3) * zmskT(3) 
     242         zw1d(4,:) = tmask(iiT+1,ijT+1,:) * zwT(4) * zmskT(4) 
     243         ptoce(:) = icb_utl_bilin_h( toce_e , iiT, ijT, zw1d ) 
     244      END IF 
     245      ! 
    231246      IF ( PRESENT(pe3t)  ) pe3t(:)  = e3t_e(iiT,ijT,:)    ! as in Nacho tarball need to be fix once we are able to reproduce Nacho results 
    232247      ! 
     
    300315         IF ( ierr > 0 ) THEN 
    301316            WRITE(numout,*) 'bottom left corner T point out of bound' 
    302             WRITE(numout,*) kii, mig( 1 ), mig(jpi) 
    303             WRITE(numout,*) kij, mjg( 1 ), mjg(jpj) 
     317            WRITE(numout,*) pi, kii, mig( 1 ), mig(jpi) 
     318            WRITE(numout,*) pj, kij, mjg( 1 ), mjg(jpj) 
     319            WRITE(numout,*) pmsk 
    304320            CALL ctl_stop('STOP','icb_utl_bilin_h: an icebergs coordinates is out of valid range (out of bound error)') 
    305321         END IF 
     
    381397      !!---------------------------------------------------------------------- 
    382398      REAL(wp), DIMENSION(0:jpi+1,0:jpj+1, jpk), INTENT(in) ::   pfld      ! field to be interpolated 
    383       REAL(wp), DIMENSION(4)                   , INTENT(in) ::   pw        ! weight 
     399      REAL(wp), DIMENSION(4,jpk)               , INTENT(in) ::   pw        ! weight 
    384400      INTEGER ,                                  INTENT(in) ::   pii, pij  ! bottom left corner 
    385401      REAL(wp), DIMENSION(jpk) :: icb_utl_bilin_3d_h 
    386402      ! 
    387403      REAL(wp), DIMENSION(4,jpk) :: zdat ! input data 
     404      INTEGER :: jk 
    388405      !!---------------------------------------------------------------------- 
    389406      ! 
     
    395412      ! 
    396413      ! compute interpolated value 
    397       icb_utl_bilin_3d_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))  
     414      DO jk=1,jpk 
     415         icb_utl_bilin_3d_h(jk) =   ( zdat(1,jk)*pw(1,jk) + zdat(2,jk)*pw(2,jk) + zdat(3,jk)*pw(3,jk) + zdat(4,jk)*pw(4,jk) ) & 
     416            &                     /   MAX(1.e-20, pw(1,jk)+pw(2,jk)+pw(3,jk)+pw(4,jk))  
     417      END DO 
    398418      ! 
    399419   END FUNCTION icb_utl_bilin_3d_h 
Note: See TracChangeset for help on using the changeset viewer.