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 15088 – NEMO

Changeset 15088


Ignore:
Timestamp:
2021-07-06T15:03:34+02:00 (3 years ago)
Author:
acc
Message:

Check in ICB changes from 2021/ticket2696_icb_halo1_halo2_compatibility (#2696) to restore halo1 - halo2 equivalence with ORCA2_ICE_PISCES and icebergs. The benefit won't be apparent until the reproducibility issue with icedyn_rhg_evp.F90 is sorted. This fix closes #2696

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

Legend:

Unmodified
Added
Removed

  • NEMO/trunk/src/OCE/ICB/icbclv.F90

    r14030 r15088  
    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/trunk/src/OCE/ICB/icbdyn.F90

    r14400 r15088  
    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/trunk/src/OCE/ICB/icbini.F90

    r14433 r15088  
    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/trunk/src/OCE/ICB/icblbc.F90

    r14433 r15088  
    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/trunk/src/OCE/ICB/icbrst.F90

    r13286 r15088  
    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/trunk/src/OCE/ICB/icbthm.F90

    r14773 r15088  
    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/trunk/src/OCE/ICB/icbutl.F90

    r14400 r15088  
    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.