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 13540 for NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/traqsr.F90 – NEMO

Ignore:
Timestamp:
2020-09-29T12:41:06+02:00 (3 years ago)
Author:
andmirek
Message:

Ticket #2386: update to latest trunk

Location:
NEMO/branches/2020/r12377_ticket2386
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/r12377_ticket2386

    • Property svn:externals
      •  

        old new  
        33^/utils/build/mk@HEAD         mk 
        44^/utils/tools@HEAD            tools 
        5 ^/vendors/AGRIF/dev@HEAD      ext/AGRIF 
         5^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 
        88 
        99# SETTE 
        10 ^/utils/CI/sette@HEAD         sette 
         10^/utils/CI/sette@13507        sette 
  • NEMO/branches/2020/r12377_ticket2386/src/OCE/TRA/traqsr.F90

    r12511 r13540  
    6363   REAL(wp) ::   xsi1r   ! inverse of rn_si1 
    6464   ! 
    65    REAL(wp) , DIMENSION(3,61)           ::   rkrgb    ! tabulated attenuation coefficients for RGB absorption 
     65   REAL(wp) , PUBLIC, DIMENSION(3,61)   ::   rkrgb    ! tabulated attenuation coefficients for RGB absorption 
    6666   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read) 
    6767 
    6868   !! * Substitutions 
    6969#  include "do_loop_substitute.h90" 
     70#  include "domzgr_substitute.h90" 
    7071   !!---------------------------------------------------------------------- 
    7172   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    110111      REAL(wp) ::   zc0 , zc1 , zc2 , zc3    !    -         - 
    111112      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 
     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 
    118118      !!---------------------------------------------------------------------- 
    119119      ! 
     
    138138            IF(lwp) WRITE(numout,*) '          nit000-1 qsr tracer content forcing field read in the restart file' 
    139139            z1_2 = 0.5_wp 
    140             CALL iom_get( numror, jpdom_autoglo, 'qsr_hc_b', qsr_hc_b, ldxios = lrxios )   ! before heat content trend due to Qsr flux 
     140            CALL iom_get( numror, jpdom_auto, 'qsr_hc_b', qsr_hc_b, ldxios = lrxios )   ! before heat content trend due to Qsr flux 
    141141         ELSE                                           ! No restart or restart not found: Euler forward time stepping 
    142142            z1_2 = 1._wp 
     
    160160      CASE( np_RGB , np_RGBc )         !==  R-G-B fluxes  ==! 
    161161         ! 
    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)   )  
     162         ALLOCATE( ze0 (jpi,jpj)           , ze1 (jpi,jpj) ,  & 
     163            &      ze2 (jpi,jpj)           , ze3 (jpi,jpj) ,  & 
     164            &      ztmp3d(jpi,jpj,nksr + 1)                     ) 
    165165         ! 
    166166         IF( nqsr == np_RGBc ) THEN          !*  Variable Chlorophyll 
    167167            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( 0, 0, 0, 0 ) 
     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( 0, 0, 0, 0, 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 
    168217            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 
     218               ztmp3d(:,:,jk) = zlui  
    190219            END DO 
    191          ELSE                                !* constant chrlorophyll 
    192            DO jk = 1, nksr + 1 
    193               zchl3d(:,:,jk) = 0.05  
    194             ENDDO 
    195220         ENDIF 
    196221         ! 
    197222         zcoef  = ( 1. - rn_abs ) / 3._wp    !* surface equi-partition in R-G-B 
    198          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         DO_2D( 0, 0, 0, 0 ) 
     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) 
    204231         END_2D 
    205232         ! 
    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 
    227          ! 
    228          DO_3D_00_00( 1, nksr ) 
    229             qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) ) 
     233         !                                    !* interior equi-partition in R-G-B depending on vertical profile of Chl 
     234         DO_3D( 0, 0, 0, 0, 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) 
    230246         END_3D 
    231247         ! 
    232          DEALLOCATE( zekb , zekg , zekr , ze0 , ze1 , ze2 , ze3 , zea , zchl3d )  
     248         DO_3D( 0, 0, 0, 0, 1, nksr )          !* now qsr induced heat content 
     249            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( ztmp3d(ji,jj,jk) - ztmp3d(ji,jj,jk+1) ) 
     250         END_3D 
     251         ! 
     252         DEALLOCATE( ze0 , ze1 , ze2 , ze3 , ztmp3d )  
    233253         ! 
    234254      CASE( np_2BD  )            !==  2-bands fluxes  ==! 
     
    236256         zz0 =        rn_abs   * r1_rho0_rcp      ! surface equi-partition in 2-bands 
    237257         zz1 = ( 1. - rn_abs ) * r1_rho0_rcp 
    238          DO_3D_00_00( 1, nksr ) 
     258         DO_3D( 0, 0, 0, 0, 1, nksr )             ! solar heat absorbed at T-point in the top 400m  
    239259            zc0 = zz0 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi1r ) 
    240260            zc1 = zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi1r ) 
     
    245265      ! 
    246266      !                          !-----------------------------! 
    247       DO_3D_00_00( 1, nksr ) 
     267      !                          !  update to the temp. trend  ! 
     268      !                          !-----------------------------! 
     269      DO_3D( 0, 0, 0, 0, 1, nksr ) 
    248270         pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   & 
    249             &                      + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t(ji,jj,jk,Kmm) 
     271            &                      + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) )   & 
     272            &                             / e3t(ji,jj,jk,Kmm) 
    250273      END_3D 
    251274      ! 
    252275      ! sea-ice: store the 1st ocean level attenuation coefficient 
    253       DO_2D_00_00 
     276      DO_2D( 0, 0, 0, 0 ) 
    254277         IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) ) 
    255278         ELSE                             ;   fraqsr_1lev(ji,jj) = 1._wp 
     
    396419         IF( .NOT.lk_top )   CALL ctl_stop( 'No bio model : ln_qsr_bio = true impossible ' ) 
    397420         ! 
     421         CALL trc_oce_rgb( rkrgb )                 ! tabulated attenuation coef. 
     422         !                                    
     423         nksr = trc_oce_ext_lev( r_si2, 33._wp )   ! level of light extinction 
     424         ! 
     425         IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m' 
     426         ! 
    398427      END SELECT 
    399428      ! 
     
    402431      ! 1st ocean level attenuation coefficient (used in sbcssm) 
    403432      IF( iom_varid( numror, 'fraqsr_1lev', ldstop = .FALSE. ) > 0 ) THEN 
    404          CALL iom_get( numror, jpdom_autoglo, 'fraqsr_1lev'  , fraqsr_1lev, ldxios = lrxios  ) 
     433         CALL iom_get( numror, jpdom_auto, 'fraqsr_1lev'  , fraqsr_1lev, ldxios = lrxios  ) 
    405434      ELSE 
    406435         fraqsr_1lev(:,:) = 1._wp   ! default : no penetration 
Note: See TracChangeset for help on using the changeset viewer.