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 13217 for NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA – NEMO

Ignore:
Timestamp:
2020-07-02T11:27:33+02:00 (4 years ago)
Author:
smasson
Message:

better e3: update with trunk@13215 see #2385

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traqsr.F90

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