Changeset 15043 for NEMO/branches/2021/ticket2696_icb_halo1_halo2_compatibility/src/OCE/ICB/icbdyn.F90

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.

File:

• NEMO/branches/2021/ticket2696_icb_halo1_halo2_compatibility/src/OCE/ICB/icbdyn.F90 (modified) (4 diffs)

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

@@ -192 +192 @@
       ld_bounced = .FALSE.
       !
-      ii0 = INT( pi0+0.5 )   ;   ij0 = INT( pj0+0.5 )       ! initial gridpoint position (T-cell)
-      ii  = INT( pi +0.5 )   ;   ij  = INT( pj +0.5 )       ! current     -         -
+      ii0 = INT( pi0+0.5 ) + (nn_hls-1)   ;   ij0 = INT( pj0+0.5 ) + (nn_hls-1)      ! initial gridpoint position (T-cell)
+      ii  = INT( pi +0.5 ) + (nn_hls-1)   ;   ij  = INT( pj +0.5 ) + (nn_hls-1)      ! current     -         -
       !
       IF( ii == ii0  .AND.  ij == ij0  )   RETURN           ! berg remains in the same cell
@@ -314 +314 @@
       zwmod  = zuwave*zuwave + zvwave*zvwave          ! The wave amplitude and length depend on the  current;
       !                                               ! wind speed relative to the ocean. Actually wmod is wmod**2 here.
-      zampl        = 0.5 * 0.02025 * zwmod            ! This is "a", the wave amplitude
-      zLwavelength =       0.32    * zwmod            ! Surface wave length fitted to data in table at
+      zampl        = 0.5_wp * 0.02025_wp * zwmod      ! This is "a", the wave amplitude
+      zLwavelength =       0.32_wp    * zwmod         ! Surface wave length fitted to data in table at
       !                                               ! http://www4.ncsu.edu/eos/users/c/ceknowle/public/chapter10/part2.html
-      zLcutoff     = 0.125 * zLwavelength
-      zLtop        = 0.25  * zLwavelength
-      zCr          = pp_Cr0 * MIN(  MAX( 0., (zL-zLcutoff) / ((zLtop-zLcutoff)+1.e-30)) , 1.)  ! Wave radiation coefficient
+      zLcutoff     = 0.125_wp * zLwavelength
+      zLtop        = 0.25_wp  * zLwavelength
+      zCr          = pp_Cr0 * MIN(  MAX( 0._wp, (zL-zLcutoff) / ((zLtop-zLcutoff)+1.e-30)) , 1._wp)  ! Wave radiation coefficient
       !                                               ! fitted to graph from Carrieres et al.,  POAC Drift Model.
-      zwave_rad    = 0.5 * pp_rho_seawater / zM * zCr * grav * zampl * MIN( zampl,zF ) * (2.*zW*zL) / (zW+zL)
+      zwave_rad    = 0.5_wp * pp_rho_seawater / zM * zCr * grav * zampl * MIN( zampl,zF ) * (2._wp*zW*zL) / (zW+zL)
       zwmod        = SQRT( zua*zua + zva*zva )        ! Wind speed
       IF( zwmod /= 0._wp ) THEN
@@ -327 +327 @@
          zvwave = zva/zwmod
       ELSE
-         zuwave = 0.   ;    zvwave=0.   ;    zwave_rad=0. ! ... and only when wind is present.     !!gm  wave_rad=0. is useless
+         zuwave = 0._wp   ;    zvwave=0._wp   ;    zwave_rad=0._wp ! ... and only when wind is present.     !!gm  wave_rad=0. is useless
       ENDIF
 
       ! Weighted drag coefficients
-      z_ocn = pp_rho_seawater / zM * (0.5*pp_Cd_wv*zW*(zD_hi)+pp_Cd_wh*zW*zL)
-      z_atm = pp_rho_air      / zM * (0.5*pp_Cd_av*zW*zF     +pp_Cd_ah*zW*zL)
-      z_ice = pp_rho_ice      / zM * (0.5*pp_Cd_iv*zW*zhi              )
+      z_ocn = pp_rho_seawater / zM * (0.5_wp*pp_Cd_wv*zW*(zD_hi)+pp_Cd_wh*zW*zL)
+      z_atm = pp_rho_air      / zM * (0.5_wp*pp_Cd_av*zW*zF     +pp_Cd_ah*zW*zL)
+      z_ice = pp_rho_ice      / zM * (0.5_wp*pp_Cd_iv*zW*zhi              )
       IF( abs(zui) + abs(zvi) == 0._wp )   z_ice = 0._wp
 
@@ -358 +358 @@
       DO itloop = 1, 2  ! Iterate on drag coefficients
          !
-         zus = 0.5 * ( zuveln + puvel )
-         zvs = 0.5 * ( zvveln + pvvel )
+         zus = 0.5_wp * ( zuveln + puvel )
+         zvs = 0.5_wp * ( zvveln + pvvel )
          zdrag_ocn = z_ocn * SQRT( (zus-zuo)*(zus-zuo) + (zvs-zvo)*(zvs-zvo) )
          zdrag_atm = z_atm * SQRT( (zus-zua)*(zus-zua) + (zvs-zva)*(zvs-zva) )