Changes between Version 2 and Version 3 of 2020WP/ENHANCE-10_acc_fix_traqsr

Timestamp:

2020-05-13T19:22:26+02:00 (4 years ago)

Author:

acc

Comment:

--

2020WP/ENHANCE-10_acc_fix_traqsr

@@ -38 +38 @@
          !
          ! code to set zchl3d(:,:,1:nskr+1)
-         !
+         ! (the chlorophyll value projected from the surface values)
          !
          zcoef  = ( 1. - rn_abs ) / 3._wp    !* surface equi-partition in R-G-B
@@ -81 +81 @@
 
 === Option 1: Minmum memory usage
-By rearranging the loop order and placing the vertical loop innermost then the code can be greatly simplified to an equivalent using minimal temporary storage:
+By rearranging the loop order and placing the vertical loop innermost then the code can be greatly simplified to an equivalent using minimal temporary storage. A couple of other changes seem sensible too, namely:
+ * rename the zchl3d array to ztmp3d (since it is now used for two purposes)
+ * only allocate ztmp3d to nksr+1; values below this are not used and nksr + 1 is likely << jpk
+ * calculate and store the attenuation coefficient look-up table index as soon as the sub-surface chlorophyll value is known. This keeps all LOG operations in one loop.
 
 {{{
       CASE( np_RGB , np_RGBc )         !==  R-G-B fluxes  ==!
          !
-         ALLOCATE( zchl3d(jpi,jpj,jpk)   )
-         !
-         ! code to set zchl3d(:,:,1:nskr+1)
-         !
-         !
+         ALLOCATE( ztmp3d(jpi,jpj,nksr + 1)   )
+         !
+         ! code to set ztmp3d(:,:,1:nskr+1)
+
+         ! including the following changes after
+         IF( nqsr == np_RGBc ) THEN          !*  Variable Chlorophyll
+            .
+            DO_3D_00_00( 1, nksr + 1 )
+            .
+            . no change until after
+               zCze    = 1.12  * (zchl)**0.803
+               !
+               ! NB. make sure zchl value is such that: zchl = MIN( 10. , MAX( 0.03, zchl ) )
+               zchl = MIN( 10. , MAX( 0.03, zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) / zdelpsi )**2 ) ) ) )
+               ! Convert chlorophyll value to attenuation coefficient look-up table index
+               ztmp3d(ji,jj,jk) = 41 + 20.*LOG10(zchl) + 1.e-15
+            END_3D
+         ELSE                                !* constant chrlorophyll
+           zchl = 0.05
+           ! NB. make sure constant value is such that:
+           zchl = MIN( 10. , MAX( 0.03, zchl ) )
+           ! Convert chlorophyll value to attenuation coefficient look-up table index
+           zlui = 41 + 20.*LOG10(zchl) + 1.e-15
+           DO jk = 1, nksr + 1
+              ztmp3d(:,:,jk) = zlui
+           END DO
+         ENDIF
+         ! 
          !
          zcoef  = ( 1. - rn_abs ) / 3._wp    !* surface equi-partition in R-G-B
-         ! store the surface SW radiation;
-         ! re-use the surface zchl3d array since the surface chl value is not used
-         zchl3d(:,:,1) = qsr(:,:)
+         ! store the surface SW radiation; re-purpose the surface ztmp3d array
+         ! since the surface attenuation coefficient is not used
+         ztmp3d(:,:,1) = qsr(:,:)
          !
          !* interior equi-partition in R-G-B depending of vertical profile of Chl
@@ -105 +131 @@
             zc4 = e3t(ji,jj,1,Kmm)
             DO jk = 2, nksr+1
-               zchl = MIN( 10. , MAX( 0.03, zchl3d(ji,jj,jk) ) )
-               irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
+               irgb = NINT( ztmp3d(ji,jj,jk) )
                zc0 = zc0 * EXP( - zc4 * xsi0r       )
                zc1 = zc1 * EXP( - zc4 * rkrgb(1,irgb) )
@@ -113 +138 @@
                zc4 = e3t(ji,jj,jk,Kmm)
                ! store the SW radiation penetrating to this location
-               ! re-use the zchl3d array since the chl value at this point will not be needed again
-               zchl3d(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk)
+               ! re-purpose the ztmp3d array since the attenuation coefficient
+               ! at this point will not be needed again
+               ztmp3d(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk)
             END DO
          END_2D
          !
          DO_3D_00_00( 1, nksr )
-            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zchl3d(ji,jj,jk) - zchl3d(ji,jj,jk+1) )
-         END_3D
-         !
-         DEALLOCATE( zchl3d )
-}}}
-
-This is code and memory efficient but will perform poorly due to non-contiguous access to the array elements (see performance section below).
-
-=== Option 2: Reduce full-depth arrays to single level arrays where possible.
-A compromise solution, which reduces memory use and maintains performance is to remove all unnecessary full-depth arrays but maintain loop order. 
+            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( ztmp3d(ji,jj,jk) - ztmp3d(ji,jj,jk+1) )
+         END_3D
+         !
+         DEALLOCATE( ztmp3d )
+}}}
+
+This is code and memory efficient but may perform poorly due to non-contiguous access to the array elements (see performance section below).
+
+=== Option 2: Low memory use (retain loop order).
+A compromise solution, which reduces memory use and should perform better is to remove all unnecessary full-depth arrays but maintain loop order by keeping a few 2D arrays. 
 {{{
        CASE( np_RGB , np_RGBc )         !==  R-G-B fluxes  ==!
          !
-         ALLOCATE( zeka(jpi,jpj)       , zekb(jpi,jpj)   ,            &
-            &      zekg(jpi,jpj)       , zekr(jpi,jpj)   ,            &
-            &      ze0 (jpi,jpj)       , ze1 (jpi,jpj) ,            &
-            &      ze2 (jpi,jpj)       , ze3 (jpi,jpj) ,            &
-            &      zchl3d(jpi,jpj,jpk)   )
-         !
-         ! code to set zchl3d(:,:,1:nskr+1)
-         !
+         ALLOCATE( ze0 (jpi,jpj)           , ze1 (jpi,jpj) ,   &
+            &      ze2 (jpi,jpj)           , ze3 (jpi,jpj) ,   &
+            &      ztmp3d(jpi,jpj,nksr + 1)                     )
+         !
+         ! code to set ztmp3d(:,:,1:nskr+1)
+         ! including the following changes after
+         IF( nqsr == np_RGBc ) THEN          !*  Variable Chlorophyll
+            .
+            DO_3D_00_00( 1, nksr + 1 )
+            .
+            . no change until after
+               zCze    = 1.12  * (zchl)**0.803
+               !
+               ! NB. make sure zchl value is such that: zchl = MIN( 10. , MAX( 0.03, zchl ) )
+               zchl = MIN( 10. , MAX( 0.03, zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) / zdelpsi )**2 ) ) ) )
+               ! Convert chlorophyll value to attenuation coefficient look-up table index
+               ztmp3d(ji,jj,jk) = 41 + 20.*LOG10(zchl) + 1.e-15
+            END_3D
+         ELSE                                !* constant chlorophyll
+            zchl = 0.05
+            ! NB. make sure constant value is such that:
+            zchl = MIN( 10. , MAX( 0.03, zchl ) )
+            ! Convert chlorophyll value to attenuation coefficient look-up table index
+            zlui = 41 + 20.*LOG10(zchl) + 1.e-15
+            DO jk = 1, nksr + 1
+               ztmp3d(:,:,jk) = zlui
+            END DO
+         ENDIF
          !
          zcoef  = ( 1. - rn_abs ) / 3._wp    !* surface equi-partition in R-G-B
@@ -147 +193 @@
             ze2(ji,jj) = zcoef  * qsr(ji,jj)
             ze3(ji,jj) = zcoef  * qsr(ji,jj)
-            ! store the surface SW radiation
-            ! re-use the surface zchl3d array since the surface chl is not used
-            zchl3d(ji,jj,1) =       qsr(ji,jj)
+            ! store the surface SW radiation; re-use the surface ztmp3d array
+            ! since the surface attenuation coefficient is not used
+            ztmp3d(ji,jj,1) =       qsr(ji,jj)
          END_2D
          !
-         DO jk = 2, nksr+1                   !* interior equi-partition in R-G-B depending of vertical profile of Chl
-            DO_2D_00_00
-               zchl = MIN( 10. , MAX( 0.03, zchl3d(ji,jj,jk) ) )
-               irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
-               ze3t = e3t(ji,jj,jk-1,Kmm)
-               zeka(ji,jj) = EXP( - ze3t * xsi0r )
-               zekb(ji,jj) = EXP( - ze3t * rkrgb(1,irgb) )
-               zekg(ji,jj) = EXP( - ze3t * rkrgb(2,irgb) )
-               zekr(ji,jj) = EXP( - ze3t * rkrgb(3,irgb) )
-            END_2D
-
-            DO_2D_00_00
-               ze0(ji,jj) = ze0(ji,jj) * zeka(ji,jj)
-               ze1(ji,jj) = ze1(ji,jj) * zekb(ji,jj)
-               ze2(ji,jj) = ze2(ji,jj) * zekg(ji,jj)
-               ze3(ji,jj) = ze3(ji,jj) * zekr(ji,jj)
-               ! store the SW radiation penetrating to this location
-               ! re-use the zchl3d array since the chl value at this point will
-               ! not be needed again
-               zchl3d(ji,jj,jk) = ( ze0(ji,jj) + ze1(ji,jj) + ze2(ji,jj) + ze3(ji,jj) ) * wmask(ji,jj,jk)
-            END_2D
-         END DO
+          !* interior equi-partition in R-G-B depending of vertical profile of Chl
+         DO_3D_00_00( 2, nksr+1 )
+            irgb = NINT( ztmp3d(ji,jj,jk) )
+            ze3t = e3t(ji,jj,jk-1,Kmm)
+            ze0(ji,jj) = ze0(ji,jj) * EXP( - ze3t * xsi0r )
+            ze1(ji,jj) = ze1(ji,jj) * EXP( - ze3t * rkrgb(1,irgb) )
+            ze2(ji,jj) = ze2(ji,jj) * EXP( - ze3t * rkrgb(2,irgb) )
+            ze3(ji,jj) = ze3(ji,jj) * EXP( - ze3t * rkrgb(3,irgb) )
+            ! store the SW radiation penetrating to this location
+            ! re-use the ztmp3d array since the attenuation coefficient
+            ! at this point will not be needed again
+            ztmp3d(ji,jj,jk) = ( ze0(ji,jj) + ze1(ji,jj) + ze2(ji,jj) + ze3(ji,jj) ) * wmask(ji,jj,jk)
+         END_3D
          !
          DO_3D_00_00( 1, nksr )
-            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zchl3d(ji,jj,jk) - zchl3d(ji,jj,jk+1) )
-         END_3D
-         !
-         DEALLOCATE( zeka, zekb , zekg , zekr , ze0 , ze1 , ze2 , ze3 , zchl3d )
-}}}
+            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( ztmp3d(ji,jj,jk) - ztmp3d(ji,jj,jk+1) )
+         END_3D
+         !
+         DEALLOCATE( ze0 , ze1 , ze2 , ze3 , ztmp3d )
+}}}
+
+=== Performance and evaluation
+Both these options produce identical results to the original code (based on an ORCA2_ICE_PISCES test using SETTE (which includes variable surface chlorophyll inputs. ln_timing was activated and the 
 ''...''