Changeset 12956


Ignore:
Timestamp:
2020-05-20T16:24:20+02:00 (4 months ago)
Author:
acc
Message:

Branch: 2020/dev_r12953_ENHANCE-10_acc_fix_traqsr. Optimised version of traqsr.F90 as discussed in wiki:2020WP/ENHANCE-10_acc_fix_traqsr (#2469). This version uses less memory than the original and achieves the same results faster. Successfully SETTE tested.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r12953_ENHANCE-10_acc_fix_traqsr/src/OCE/TRA/traqsr.F90

    r12489 r12956  
    110110      REAL(wp) ::   zc0 , zc1 , zc2 , zc3    !    -         - 
    111111      REAL(wp) ::   zzc0, zzc1, zzc2, zzc3   !    -         - 
    112       REAL(wp) ::   zz0 , zz1                !    -         - 
    113       REAL(wp) ::   zCb, zCmax, zze, zpsi, zpsimax, zdelpsi, zCtot, zCze 
    114       REAL(wp) ::   zlogc, zlogc2, zlogc3  
    115       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   :: zekb, zekg, zekr 
    116       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt 
    117       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zetot, zchl3d 
     112      REAL(wp) ::   zz0 , zz1 , ze3t, zlui   !    -         - 
     113      REAL(wp) ::   zCb, zCmax, zpsi, zpsimax, zrdpsi, zCze 
     114      REAL(wp) ::   zlogc, zlogze, zlogCtot, zlogCze 
     115      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   :: ze0, ze1, ze2, ze3 
     116      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt, zetot, ztmp3d 
    118117      !!---------------------------------------------------------------------- 
    119118      ! 
     
    160159      CASE( np_RGB , np_RGBc )         !==  R-G-B fluxes  ==! 
    161160         ! 
    162          ALLOCATE( zekb(jpi,jpj)     , zekg(jpi,jpj)     , zekr  (jpi,jpj)     , & 
    163             &      ze0 (jpi,jpj,jpk) , ze1 (jpi,jpj,jpk) , ze2   (jpi,jpj,jpk) , & 
    164             &      ze3 (jpi,jpj,jpk) , zea (jpi,jpj,jpk) , zchl3d(jpi,jpj,jpk)   )  
     161         ALLOCATE( ze0 (jpi,jpj)           , ze1 (jpi,jpj) ,  & 
     162            &      ze2 (jpi,jpj)           , ze3 (jpi,jpj) ,  & 
     163            &      ztmp3d(jpi,jpj,nksr + 1)                     ) 
    165164         ! 
    166165         IF( nqsr == np_RGBc ) THEN          !*  Variable Chlorophyll 
    167166            CALL fld_read( kt, 1, sf_chl )         ! Read Chl data and provides it at the current time step 
     167            ! 
     168            ! Separation in R-G-B depending on the surface Chl 
     169            ! perform and store as many of the 2D calculations as possible 
     170            ! before the 3D loop (use the temporary 2D arrays to replace the 
     171            ! most expensive calculations) 
     172            ! 
     173            ze0(:,:) = LOG ( MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(:,:,1) ) ) )     ! ze0 = log(zchl) 
     174            ze1(:,:) = 0.113328685307 + 0.803 * ze0(:,:)                           ! ze1 : log(zCze)  = log (1.12  * zchl**0.803) 
     175            ze2(:,:) = 3.703768066608 + 0.459 * ze0(:,:)                           ! ze2 : log(zCtot) = log(40.6  * zchl**0.459) 
     176            ze3(:,:) = 6.34247346942  - 0.746 * ze2(:,:)                           ! ze3 : log(zze)   = log(568.2 * zCtot**(-0.746)) 
     177            WHERE( ze3 > 4.62497281328 ) ze3 = 5.298317366548 - 0.293 * ze2        !       IF( log(zze) > log(102.) ) log(zze) =  
     178            !                                                                      !                    log(200.0 * zCtot**(-0.293)) 
     179            ze1(:,:) = EXP(ze1(:,:))                                               ! ze1 = zCze 
     180            ze2(:,:) = 1._wp / ( 0.710 + ze0(:,:) * ( 0.159 + ze0(:,:) * 0.021 ) ) ! ze2 = 1/zdelpsi 
     181            ze3(:,:) = 1._wp / EXP(ze3(:,:))                                       ! ze3 = 1/zze 
     182            ! 
     183            DO_3D_00_00 ( 1, nksr + 1 ) 
     184               ! zchl    = ALOG( ze0(ji,jj) ) 
     185               zlogc = ze0(ji,jj) 
     186               ! 
     187               zCb       = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) ) 
     188               zCmax     = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 ) 
     189               zpsimax   = 0.6   - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) ) 
     190               ! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 ) 
     191               ! 
     192               zCze   = ze1(ji,jj) 
     193               zrdpsi = ze2(ji,jj)                                                 ! 1/zdelpsi 
     194               zpsi   = ze3(ji,jj) * gdepw(ji,jj,jk,Kmm)                           ! gdepw/zze 
     195               ! 
     196               ! NB. make sure zchl value is such that: zchl = MIN( 10. , MAX( 0.03, zchl ) ) 
     197               zchl = MIN( 10. , MAX( 0.03, zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) ) ) 
     198               ! Convert chlorophyll value to attenuation coefficient look-up table index 
     199               ztmp3d(ji,jj,jk) = 41 + 20.*LOG10(zchl) + 1.e-15 
     200            END_3D 
     201         ELSE                                !* constant chlorophyll 
     202            zchl = 0.05 
     203            ! NB. make sure constant value is such that:  
     204            zchl = MIN( 10. , MAX( 0.03, zchl ) ) 
     205            ! Convert chlorophyll value to attenuation coefficient look-up table index 
     206            zlui = 41 + 20.*LOG10(zchl) + 1.e-15 
    168207            DO jk = 1, nksr + 1 
    169                DO jj = 2, jpjm1                       ! Separation in R-G-B depending of the surface Chl 
    170                   DO ji = 2, jpim1 
    171                      zchl    = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) 
    172                      zCtot   = 40.6  * zchl**0.459 
    173                      zze     = 568.2 * zCtot**(-0.746) 
    174                      IF( zze > 102. ) zze = 200.0 * zCtot**(-0.293) 
    175                      zpsi    = gdepw(ji,jj,jk,Kmm) / zze 
    176                      ! 
    177                      zlogc   = LOG( zchl ) 
    178                      zlogc2  = zlogc * zlogc 
    179                      zlogc3  = zlogc * zlogc * zlogc 
    180                      zCb     = 0.768 + 0.087 * zlogc - 0.179 * zlogc2 - 0.025 * zlogc3 
    181                      zCmax   = 0.299 - 0.289 * zlogc + 0.579 * zlogc2 
    182                      zpsimax = 0.6   - 0.640 * zlogc + 0.021 * zlogc2 + 0.115 * zlogc3 
    183                      zdelpsi = 0.710 + 0.159 * zlogc + 0.021 * zlogc2 
    184                      zCze    = 1.12  * (zchl)**0.803  
    185                      ! 
    186                      zchl3d(ji,jj,jk) = zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) / zdelpsi )**2 ) ) 
    187                   END DO 
    188                   ! 
    189                END DO 
     208               ztmp3d(:,:,jk) = zlui  
    190209            END DO 
    191          ELSE                                !* constant chrlorophyll 
    192            DO jk = 1, nksr + 1 
    193               zchl3d(:,:,jk) = 0.05  
    194             ENDDO 
    195210         ENDIF 
    196211         ! 
    197212         zcoef  = ( 1. - rn_abs ) / 3._wp    !* surface equi-partition in R-G-B 
    198213         DO_2D_00_00 
    199             ze0(ji,jj,1) = rn_abs * qsr(ji,jj) 
    200             ze1(ji,jj,1) = zcoef  * qsr(ji,jj) 
    201             ze2(ji,jj,1) = zcoef  * qsr(ji,jj) 
    202             ze3(ji,jj,1) = zcoef  * qsr(ji,jj) 
    203             zea(ji,jj,1) =          qsr(ji,jj) 
     214            ze0(ji,jj) = rn_abs * qsr(ji,jj) 
     215            ze1(ji,jj) = zcoef  * qsr(ji,jj) 
     216            ze2(ji,jj) = zcoef  * qsr(ji,jj) 
     217            ze3(ji,jj) = zcoef  * qsr(ji,jj) 
     218            ! store the surface SW radiation; re-use the surface ztmp3d array  
     219            ! since the surface attenuation coefficient is not used 
     220            ztmp3d(ji,jj,1) =       qsr(ji,jj) 
    204221         END_2D 
    205222         ! 
    206          DO jk = 2, nksr+1                   !* interior equi-partition in R-G-B depending of vertical profile of Chl 
    207             DO_2D_00_00 
    208                zchl = MIN( 10. , MAX( 0.03, zchl3d(ji,jj,jk) ) ) 
    209                irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) 
    210                zekb(ji,jj) = rkrgb(1,irgb) 
    211                zekg(ji,jj) = rkrgb(2,irgb) 
    212                zekr(ji,jj) = rkrgb(3,irgb) 
    213             END_2D 
    214  
    215             DO_2D_00_00 
    216                zc0 = ze0(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * xsi0r       ) 
    217                zc1 = ze1(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * zekb(ji,jj) ) 
    218                zc2 = ze2(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * zekg(ji,jj) ) 
    219                zc3 = ze3(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * zekr(ji,jj) ) 
    220                ze0(ji,jj,jk) = zc0 
    221                ze1(ji,jj,jk) = zc1 
    222                ze2(ji,jj,jk) = zc2 
    223                ze3(ji,jj,jk) = zc3 
    224                zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk) 
    225             END_2D 
    226          END DO 
     223         !* interior equi-partition in R-G-B depending on vertical profile of Chl 
     224         DO_3D_00_00 ( 2, nksr + 1 ) 
     225            ze3t = e3t(ji,jj,jk-1,Kmm) 
     226            irgb = NINT( ztmp3d(ji,jj,jk) ) 
     227            zc0 = ze0(ji,jj) * EXP( - ze3t * xsi0r ) 
     228            zc1 = ze1(ji,jj) * EXP( - ze3t * rkrgb(1,irgb) ) 
     229            zc2 = ze2(ji,jj) * EXP( - ze3t * rkrgb(2,irgb) ) 
     230            zc3 = ze3(ji,jj) * EXP( - ze3t * rkrgb(3,irgb) ) 
     231            ze0(ji,jj) = zc0 
     232            ze1(ji,jj) = zc1 
     233            ze2(ji,jj) = zc2 
     234            ze3(ji,jj) = zc3 
     235            ztmp3d(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk) 
     236         END_3D 
    227237         ! 
    228238         DO_3D_00_00( 1, nksr ) 
    229             qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) ) 
     239            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( ztmp3d(ji,jj,jk) - ztmp3d(ji,jj,jk+1) ) 
    230240         END_3D 
    231241         ! 
    232          DEALLOCATE( zekb , zekg , zekr , ze0 , ze1 , ze2 , ze3 , zea , zchl3d )  
     242         DEALLOCATE( ze0 , ze1 , ze2 , ze3 , ztmp3d )  
    233243         ! 
    234244      CASE( np_2BD  )            !==  2-bands fluxes  ==! 
Note: See TracChangeset for help on using the changeset viewer.