Changeset 13205


Ignore:
Timestamp:
2020-07-02T10:43:08+02:00 (5 weeks ago)
Author:
acc
Message:

Mid-year merge. Merge changes from 2020/dev_r12953_ENHANCE-10_acc_fix_traqsr onto trunk (traqsr.F90 only)

Location:
NEMO/trunk
Files:
1 edited
1 copied

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/OCE/TRA/traqsr.F90

    r12489 r13205  
    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            DO_2D_00_00 
     174                       ! zlogc = log(zchl) 
     175               zlogc = LOG ( MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) )      
     176                       ! zc1 : log(zCze)  = log (1.12  * zchl**0.803) 
     177               zc1   = 0.113328685307 + 0.803 * zlogc                          
     178                       ! zc2 : log(zCtot) = log(40.6  * zchl**0.459) 
     179               zc2   = 3.703768066608 + 0.459 * zlogc                            
     180                       ! zc3 : log(zze)   = log(568.2 * zCtot**(-0.746)) 
     181               zc3   = 6.34247346942  - 0.746 * zc2                            
     182                       ! IF( log(zze) > log(102.) ) log(zze) = log(200.0 * zCtot**(-0.293)) 
     183               IF( zc3 > 4.62497281328 ) zc3 = 5.298317366548 - 0.293 * zc2         
     184               !    
     185               ze0(ji,jj) = zlogc                                                 ! ze0 = log(zchl) 
     186               ze1(ji,jj) = EXP( zc1 )                                            ! ze1 = zCze 
     187               ze2(ji,jj) = 1._wp / ( 0.710 + zlogc * ( 0.159 + zlogc * 0.021 ) ) ! ze2 = 1/zdelpsi 
     188               ze3(ji,jj) = EXP( - zc3 )                                          ! ze3 = 1/zze 
     189            END_2D 
     190             
     191! 
     192            DO_3D_00_00 ( 1, nksr + 1 ) 
     193               ! zchl    = ALOG( ze0(ji,jj) ) 
     194               zlogc = ze0(ji,jj) 
     195               ! 
     196               zCb       = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) ) 
     197               zCmax     = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 ) 
     198               zpsimax   = 0.6   - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) ) 
     199               ! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 ) 
     200               ! 
     201               zCze   = ze1(ji,jj) 
     202               zrdpsi = ze2(ji,jj)                                                 ! 1/zdelpsi 
     203               zpsi   = ze3(ji,jj) * gdepw(ji,jj,jk,Kmm)                           ! gdepw/zze 
     204               ! 
     205               ! NB. make sure zchl value is such that: zchl = MIN( 10. , MAX( 0.03, zchl ) ) 
     206               zchl = MIN( 10. , MAX( 0.03, zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) ) ) 
     207               ! Convert chlorophyll value to attenuation coefficient look-up table index 
     208               ztmp3d(ji,jj,jk) = 41 + 20.*LOG10(zchl) + 1.e-15 
     209            END_3D 
     210         ELSE                                !* constant chlorophyll 
     211            zchl = 0.05 
     212            ! NB. make sure constant value is such that:  
     213            zchl = MIN( 10. , MAX( 0.03, zchl ) ) 
     214            ! Convert chlorophyll value to attenuation coefficient look-up table index 
     215            zlui = 41 + 20.*LOG10(zchl) + 1.e-15 
    168216            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 
     217               ztmp3d(:,:,jk) = zlui  
    190218            END DO 
    191          ELSE                                !* constant chrlorophyll 
    192            DO jk = 1, nksr + 1 
    193               zchl3d(:,:,jk) = 0.05  
    194             ENDDO 
    195219         ENDIF 
    196220         ! 
    197221         zcoef  = ( 1. - rn_abs ) / 3._wp    !* surface equi-partition in R-G-B 
    198222         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) 
     223            ze0(ji,jj) = rn_abs * qsr(ji,jj) 
     224            ze1(ji,jj) = zcoef  * qsr(ji,jj) 
     225            ze2(ji,jj) = zcoef  * qsr(ji,jj) 
     226            ze3(ji,jj) = zcoef  * qsr(ji,jj) 
     227            ! store the surface SW radiation; re-use the surface ztmp3d array  
     228            ! since the surface attenuation coefficient is not used 
     229            ztmp3d(ji,jj,1) =       qsr(ji,jj) 
    204230         END_2D 
    205231         ! 
    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 
     232         !* interior equi-partition in R-G-B depending on vertical profile of Chl 
     233         DO_3D_00_00 ( 2, nksr + 1 ) 
     234            ze3t = e3t(ji,jj,jk-1,Kmm) 
     235            irgb = NINT( ztmp3d(ji,jj,jk) ) 
     236            zc0 = ze0(ji,jj) * EXP( - ze3t * xsi0r ) 
     237            zc1 = ze1(ji,jj) * EXP( - ze3t * rkrgb(1,irgb) ) 
     238            zc2 = ze2(ji,jj) * EXP( - ze3t * rkrgb(2,irgb) ) 
     239            zc3 = ze3(ji,jj) * EXP( - ze3t * rkrgb(3,irgb) ) 
     240            ze0(ji,jj) = zc0 
     241            ze1(ji,jj) = zc1 
     242            ze2(ji,jj) = zc2 
     243            ze3(ji,jj) = zc3 
     244            ztmp3d(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk) 
     245         END_3D 
    227246         ! 
    228247         DO_3D_00_00( 1, nksr ) 
    229             qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) ) 
     248            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( ztmp3d(ji,jj,jk) - ztmp3d(ji,jj,jk+1) ) 
    230249         END_3D 
    231250         ! 
    232          DEALLOCATE( zekb , zekg , zekr , ze0 , ze1 , ze2 , ze3 , zea , zchl3d )  
     251         DEALLOCATE( ze0 , ze1 , ze2 , ze3 , ztmp3d )  
    233252         ! 
    234253      CASE( np_2BD  )            !==  2-bands fluxes  ==! 
Note: See TracChangeset for help on using the changeset viewer.