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 15043 for NEMO – NEMO

Changeset 15043 for NEMO


Ignore:
Timestamp:
2021-06-22T13:26:16+02:00 (3 years ago)
Author:
acc
Message:

#2696 .A possible fix for obtaining identical results with icebergs active and with either
nn_hls=1 or nn_hls=2. This solution works by ensuring the real-valued iceberg position
calculations are done with identical numbers in both cases. Iceberg positions are stored
as a continuous local grid-cell value and a factor of (nn_hls-1) is subtracted from
initial positions to ensure equivalence in both cases. This factor needs to be added
back whenever iceberg positions are converted to local indices. The changes in this
commit are enough to ensure a full ORCA2_ICE_PISCES SETTE test is passed in both case
(992 timesteps) and provides identical results. Suitability of this solution for inclusion
in the trunk needs to be discussed and verified on other platforms and compilers.

Location:
NEMO/branches/2021/ticket2696_icb_halo1_halo2_compatibility/src/OCE/ICB
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2021/ticket2696_icb_halo1_halo2_compatibility/src/OCE/ICB/icbclv.F90

    r14030 r15043  
    133133                  newpt%lon = glamt(ji,jj)         ! at t-point (centre of the cell) 
    134134                  newpt%lat = gphit(ji,jj) 
    135                   newpt%xi  = REAL( mig(ji), wp ) 
    136                   newpt%yj  = REAL( mjg(jj), wp ) 
     135                  newpt%xi  = REAL( mig(ji), wp ) - ( nn_hls - 1 ) 
     136                  newpt%yj  = REAL( mjg(jj), wp ) - ( nn_hls - 1 ) 
    137137                  ! 
    138138                  newpt%uvel = 0._wp               ! initially at rest 

  • NEMO/branches/2021/ticket2696_icb_halo1_halo2_compatibility/src/OCE/ICB/icbdyn.F90

    r14400 r15043  
    192192      ld_bounced = .FALSE. 
    193193      ! 
    194       ii0 = INT( pi0+0.5 )   ;   ij0 = INT( pj0+0.5 )       ! initial gridpoint position (T-cell) 
    195       ii  = INT( pi +0.5 )   ;   ij  = INT( pj +0.5 )       ! current     -         - 
     194      ii0 = INT( pi0+0.5 ) + (nn_hls-1)   ;   ij0 = INT( pj0+0.5 ) + (nn_hls-1)      ! initial gridpoint position (T-cell) 
     195      ii  = INT( pi +0.5 ) + (nn_hls-1)   ;   ij  = INT( pj +0.5 ) + (nn_hls-1)      ! current     -         - 
    196196      ! 
    197197      IF( ii == ii0  .AND.  ij == ij0  )   RETURN           ! berg remains in the same cell 
    … …  
    314314      zwmod  = zuwave*zuwave + zvwave*zvwave          ! The wave amplitude and length depend on the  current; 
    315315      !                                               ! wind speed relative to the ocean. Actually wmod is wmod**2 here. 
    316       zampl        = 0.5 * 0.02025 * zwmod            ! This is "a", the wave amplitude 
    317       zLwavelength =       0.32    * zwmod            ! Surface wave length fitted to data in table at 
     316      zampl        = 0.5_wp * 0.02025_wp * zwmod      ! This is "a", the wave amplitude 
     317      zLwavelength =       0.32_wp    * zwmod         ! Surface wave length fitted to data in table at 
    318318      !                                               ! http://www4.ncsu.edu/eos/users/c/ceknowle/public/chapter10/part2.html 
    319       zLcutoff     = 0.125 * zLwavelength 
    320       zLtop        = 0.25  * zLwavelength 
    321       zCr          = pp_Cr0 * MIN(  MAX( 0., (zL-zLcutoff) / ((zLtop-zLcutoff)+1.e-30)) , 1.)  ! Wave radiation coefficient 
     319      zLcutoff     = 0.125_wp * zLwavelength 
     320      zLtop        = 0.25_wp  * zLwavelength 
     321      zCr          = pp_Cr0 * MIN(  MAX( 0._wp, (zL-zLcutoff) / ((zLtop-zLcutoff)+1.e-30)) , 1._wp)  ! Wave radiation coefficient 
    322322      !                                               ! fitted to graph from Carrieres et al.,  POAC Drift Model. 
    323       zwave_rad    = 0.5 * pp_rho_seawater / zM * zCr * grav * zampl * MIN( zampl,zF ) * (2.*zW*zL) / (zW+zL) 
     323      zwave_rad    = 0.5_wp * pp_rho_seawater / zM * zCr * grav * zampl * MIN( zampl,zF ) * (2._wp*zW*zL) / (zW+zL) 
    324324      zwmod        = SQRT( zua*zua + zva*zva )        ! Wind speed 
    325325      IF( zwmod /= 0._wp ) THEN 
    … …  
    327327         zvwave = zva/zwmod 
    328328      ELSE 
    329          zuwave = 0.   ;    zvwave=0.   ;    zwave_rad=0. ! ... and only when wind is present.     !!gm  wave_rad=0. is useless 
     329         zuwave = 0._wp   ;    zvwave=0._wp   ;    zwave_rad=0._wp ! ... and only when wind is present.     !!gm  wave_rad=0. is useless 
    330330      ENDIF 
    331331 
    332332      ! Weighted drag coefficients 
    333       z_ocn = pp_rho_seawater / zM * (0.5*pp_Cd_wv*zW*(zD_hi)+pp_Cd_wh*zW*zL) 
    334       z_atm = pp_rho_air      / zM * (0.5*pp_Cd_av*zW*zF     +pp_Cd_ah*zW*zL) 
    335       z_ice = pp_rho_ice      / zM * (0.5*pp_Cd_iv*zW*zhi              ) 
     333      z_ocn = pp_rho_seawater / zM * (0.5_wp*pp_Cd_wv*zW*(zD_hi)+pp_Cd_wh*zW*zL) 
     334      z_atm = pp_rho_air      / zM * (0.5_wp*pp_Cd_av*zW*zF     +pp_Cd_ah*zW*zL) 
     335      z_ice = pp_rho_ice      / zM * (0.5_wp*pp_Cd_iv*zW*zhi              ) 
    336336      IF( abs(zui) + abs(zvi) == 0._wp )   z_ice = 0._wp 
    337337 
    … …  
    358358      DO itloop = 1, 2  ! Iterate on drag coefficients 
    359359         ! 
    360          zus = 0.5 * ( zuveln + puvel ) 
    361          zvs = 0.5 * ( zvveln + pvvel ) 
     360         zus = 0.5_wp * ( zuveln + puvel ) 
     361         zvs = 0.5_wp * ( zvveln + pvvel ) 
    362362         zdrag_ocn = z_ocn * SQRT( (zus-zuo)*(zus-zuo) + (zvs-zvo)*(zvs-zvo) ) 
    363363         zdrag_atm = z_atm * SQRT( (zus-zua)*(zus-zua) + (zvs-zva)*(zvs-zva) ) 

  • NEMO/branches/2021/ticket2696_icb_halo1_halo2_compatibility/src/OCE/ICB/icbini.F90

    r14433 r15043  
    182182      i3 = INT( src_calving(i1,jpj/2) ) 
    183183      jj = INT( i3/nicbpack ) 
    184       ricb_left = REAL( i3 - nicbpack*jj, wp ) 
     184      ricb_left = REAL( i3 - nicbpack*jj, wp ) - (nn_hls-1) 
    185185      i1 = MIN( nicbei+1, jpi ) 
    186186      i3 = INT( src_calving(i1,jpj/2) ) 
    187187      jj = INT( i3/nicbpack ) 
    188       ricb_right = REAL( i3 - nicbpack*jj, wp ) 
     188      ricb_right = REAL( i3 - nicbpack*jj, wp ) - (nn_hls-1) 
    189189       
    190190      ! north fold 
    … …  
    360360                rn_test_box(3) < gphit(ji,jj) .AND. gphit(ji,jj) < rn_test_box(4) ) THEN 
    361361               localberg%mass_scaling = rn_mass_scaling(iberg) 
    362                localpt%xi = REAL( mig(ji), wp ) 
    363                localpt%yj = REAL( mjg(jj), wp ) 
     362               localpt%xi = REAL( mig(ji) - (nn_hls-1), wp ) 
     363               localpt%yj = REAL( mjg(jj) - (nn_hls-1), wp ) 
    364364               CALL icb_utl_interp( localpt%xi, localpt%yj, plat=localpt%lat, plon=localpt%lon )    
    365365               localpt%mass      = rn_initial_mass     (iberg) 

  • NEMO/branches/2021/ticket2696_icb_halo1_halo2_compatibility/src/OCE/ICB/icblbc.F90

    r14433 r15043  
    229229         DO WHILE (ASSOCIATED(this)) 
    230230            pt => this%current_point 
    231             IF( ipe_E >= 0 .AND. pt%xi > REAL(mig(nicbei),wp) + 0.5_wp ) THEN 
     231            IF( ipe_E >= 0 .AND. pt%xi > REAL(mig(nicbei),wp) + 0.5_wp - (nn_hls-1) ) THEN 
    232232               tmpberg => this 
    233233               this => this%next 
    … …  
    242242               CALL icb_pack_into_buffer( tmpberg, obuffer_e, ibergs_to_send_e) 
    243243               CALL icb_utl_delete(first_berg, tmpberg) 
    244             ELSE IF( ipe_W >= 0 .AND. pt%xi < REAL(mig(nicbdi),wp) - 0.5_wp ) THEN 
     244            ELSE IF( ipe_W >= 0 .AND. pt%xi < REAL(mig(nicbdi),wp) - 0.5_wp - (nn_hls-1) ) THEN 
    245245               tmpberg => this 
    246246               this => this%next 
    … …  
    321321         DO WHILE (ASSOCIATED(this)) 
    322322            pt => this%current_point 
    323             IF( ipe_N >= 0 .AND. pt%yj > REAL(mjg(nicbej),wp) + 0.5_wp ) THEN 
     323            IF( ipe_N >= 0 .AND. pt%yj > REAL(mjg(nicbej),wp) + 0.5_wp - (nn_hls-1) ) THEN 
    324324               tmpberg => this 
    325325               this => this%next 
    … …  
    331331               CALL icb_pack_into_buffer( tmpberg, obuffer_n, ibergs_to_send_n) 
    332332               CALL icb_utl_delete(first_berg, tmpberg) 
    333             ELSE IF( ipe_S >= 0 .AND. pt%yj < REAL(mjg(nicbdj),wp) - 0.5_wp ) THEN 
     333            ELSE IF( ipe_S >= 0 .AND. pt%yj < REAL(mjg(nicbdj),wp) - 0.5_wp - (nn_hls-1) ) THEN 
    334334               tmpberg => this 
    335335               this => this%next 
    … …  
    442442         DO WHILE (ASSOCIATED(this)) 
    443443            pt => this%current_point 
    444             IF( pt%xi < REAL(mig(nicbdi),wp) - 0.5_wp .OR. & 
    445                 pt%xi > REAL(mig(nicbei),wp) + 0.5_wp .OR. & 
    446                 pt%yj < REAL(mjg(nicbdj),wp) - 0.5_wp .OR. & 
    447                 pt%yj > REAL(mjg(nicbej),wp) + 0.5_wp ) THEN 
     444            IF( pt%xi < REAL(mig(nicbdi),wp) - 0.5_wp - (nn_hls-1) .OR. & 
     445                pt%xi > REAL(mig(nicbei),wp) + 0.5_wp - (nn_hls-1) .OR. & 
     446                pt%yj < REAL(mjg(nicbdj),wp) - 0.5_wp - (nn_hls-1) .OR. & 
     447                pt%yj > REAL(mjg(nicbej),wp) + 0.5_wp - (nn_hls-1) ) THEN 
    448448               i = i + 1 
    449449               WRITE(numicb,*) 'berg lost in halo: ', this%number(:) 
    … …  
    514514               DO WHILE (ASSOCIATED(this)) 
    515515                  pt => this%current_point 
    516                   iine = INT( pt%xi + 0.5 ) 
     516                  iine = INT( pt%xi + 0.5 ) + (nn_hls-1) 
    517517                  iproc = nicbflddest(mi1(iine)) 
    518                   IF( pt%yj > REAL(mjg(nicbej),wp) + 0.5_wp ) THEN 
     518                  IF( pt%yj > REAL(mjg(nicbej),wp) + 0.5_wp - (nn_hls-1) ) THEN 
    519519                     IF( iproc == ifldproc ) THEN 
    520520                        ! 
    … …  
    592592               DO WHILE (ASSOCIATED(this)) 
    593593                  pt => this%current_point 
    594                   iine = INT( pt%xi + 0.5 ) 
    595                   ijne = INT( pt%yj + 0.5 ) 
     594                  iine = INT( pt%xi + 0.5 ) + (nn_hls-1) 
     595                  ijne = INT( pt%yj + 0.5 ) + (nn_hls-1) 
    596596                  ipts  = nicbfldpts (mi1(iine)) 
    597597                  iproc = nicbflddest(mi1(iine)) 
    598                   IF( pt%yj > REAL(mjg(nicbej),wp) + 0.5_wp ) THEN 
     598                  IF( pt%yj > REAL(mjg(nicbej),wp) + 0.5_wp - (nn_hls-1) ) THEN 
    599599                     IF( iproc == ifldproc ) THEN 
    600600                        ! 

  • NEMO/branches/2021/ticket2696_icb_halo1_halo2_compatibility/src/OCE/ICB/icbrst.F90

    r13286 r15043  
    8888            CALL iom_get( ncid, 'yj'     ,localpt%yj  , ktime=jn ) 
    8989 
    90             ii = INT( localpt%xi + 0.5 ) 
    91             ij = INT( localpt%yj + 0.5 ) 
     90            ii = INT( localpt%xi + 0.5 ) + ( nn_hls-1 ) 
     91            ij = INT( localpt%yj + 0.5 ) + ( nn_hls-1 ) 
    9292            ! Only proceed if this iceberg is on the local processor (excluding halos). 
    9393            IF ( ii >= mig(Nis0) .AND. ii <= mig(Nie0) .AND.   & 

  • NEMO/branches/2021/ticket2696_icb_halo1_halo2_compatibility/src/OCE/ICB/icbthm.F90

    r14773 r15043  
    113113         zyj  = pt%yj 
    114114         ii  = INT( zxi + 0.5 )                            ! T-cell of the berg 
    115          ii  = mi1( ii ) 
     115         ii  = mi1( ii + (nn_hls-1) ) 
    116116         ij  = INT( zyj + 0.5 )               
    117          ij  = mj1( ij ) 
     117         ij  = mj1( ij + (nn_hls-1) ) 
    118118         zVol = zT * zW * zL 
    119119 
    … …  
    203203            zLbits   = MIN( zL, zW, zT, 40._wp )                                     ! assume bergy bits are smallest dimension or 40 meters 
    204204            zAbits   = ( zMbits / rn_rho_bergs ) / zLbits                            ! Effective bottom area (assuming T=Lbits) 
    205             zMbb     = MAX( 0.58_wp*(zdvo**0.8_wp)*(zSST+2._wp) /   & 
     205            zMbb     = MAX( 0.58_wp*(zdvob**0.8_wp)*(zSST+2._wp) /   & 
    206206               &                              ( zLbits**0.2_wp ) , 0._wp ) * z1_rday ! Basal turbulent melting (for bits) 
    207207            zMbb     = rn_rho_bergs * zAbits * zMbb                                  ! in kg/s 

  • NEMO/branches/2021/ticket2696_icb_halo1_halo2_compatibility/src/OCE/ICB/icbutl.F90

    r14400 r15043  
    300300         zwj = pj - 0.5_wp - REAL(kij,wp) 
    301301      END SELECT 
     302      kii = kii + (nn_hls-1) 
     303      kij = kij + (nn_hls-1) 
    302304      ! 
    303305      ! compute weight 
    … …  
    461463 
    462464      ! conversion to local domain (no need to do a sanity check already done in icbpos) 
    463       ii = mi1(ii) 
    464       ij = mj1(ij) 
     465      ii = mi1(ii) + (nn_hls-1) 
     466      ij = mj1(ij) + (nn_hls-1) 
    465467      ! 
    466468      IF(    0.0_wp <= zi .AND. zi < 0.5_wp   ) THEN 
Note: See TracChangeset for help on using the changeset viewer.