MODULE trcopt !!====================================================================== !! *** MODULE trcopt *** !! TOP : Compute the light in the water column for RGB wavelengths !!====================================================================== !! History : 1.0 ! 2020 (T. Lovato) Initial code !!---------------------------------------------------------------------- !! trc_opt : light availability in the water column !!---------------------------------------------------------------------- USE trc ! tracer variables USE oce_trc ! tracer-ocean share variables USE iom ! I/O manager USE fldread ! time interpolation IMPLICIT NONE PRIVATE PUBLIC trc_opt ! called in spefici BGC model routines PUBLIC trc_opt_ini ! called in trcini.F90 PUBLIC trc_opt_alloc !! * Shared module variables LOGICAL :: ln_varpar ! boolean for variable PAR fraction REAL(wp) :: parlux ! Fraction of shortwave as PAR CHARACTER (len=25) :: light_loc ! Light location in the water cell ('center', 'integral') TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_par ! structure of input par INTEGER :: ntimes_par ! number of time steps in par file REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: par_varsw ! PAR fraction of shortwave REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ekb, ekg, ekr ! wavelength (Red-Green-Blue) REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xeps ! weighted diffusion coefficient INTEGER :: nksrp ! levels below which the light cannot penetrate ( depth larger than 391 m) ! TL: This array should come directly from traqsr module REAL(wp), DIMENSION(3,61) :: xkrgb ! tabulated attenuation coefficients for RGB absorption !! * Substitutions # include "do_loop_substitute.h90" # include "domzgr_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/TOP 4.0 , NEMO Consortium (2020) !! $Id: trcopt.F90 12377 2020-02-12 14:39:06Z acc $ !! Software governed by the CeCILL license (see ./LICENSE) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE trc_opt( kt, knt, Kbb, Kmm, zchl, ze1, ze2, ze3) !!--------------------------------------------------------------------- !! *** ROUTINE trc_opt *** !! !! ** Purpose : Compute the light availability in the water column !! depending on depth and chlorophyll concentration !! !! ** Method : Morel !!--------------------------------------------------------------------- INTEGER, INTENT(in) :: kt, knt ! ocean time step INTEGER, INTENT(in) :: Kbb, Kmm ! time level indices REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: zchl ! chlorophyll field REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(out),OPTIONAL :: ze1, ze2, ze3 ! PAR for individual wavelength ! INTEGER :: ji, jj, jk, irgb REAL(wp) :: ztmp REAL(wp), DIMENSION(jpi,jpj ) :: parsw, zqsr100, zqsr_corr REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze0 !!--------------------------------------------------------------------- ! IF( ln_timing ) CALL timing_start('trc_opt') ! Initialisation of variables used to compute PAR ! ----------------------------------------------- ze0(:,:,:) = 0._wp ze1(:,:,:) = 0._wp ze2(:,:,:) = 0._wp ze3(:,:,:) = 0._wp ! PAR conversion factor ! -------------------- IF( knt == 1 .AND. ln_varpar ) CALL trc_opt_sbc( kt ) ! IF( ln_varpar ) THEN ; parsw(:,:) = par_varsw(:,:) ELSE ; parsw(:,:) = parlux / 3.0 ENDIF ! Attenuation coef. function of Chlorophyll and wavelength (RGB) ! -------------------------------------------------------------- DO_3D( 1, 1, 1, 1, 1, jpkm1 ) ztmp = ( zchl(ji,jj,jk) + rtrn ) * 1.e6 ztmp = MIN( 10. , MAX( 0.05, ztmp ) ) irgb = NINT( 41 + 20.* LOG10( ztmp ) + rtrn ) ! ekb(ji,jj,jk) = xkrgb(1,irgb) * e3t(ji,jj,jk,Kmm) ekg(ji,jj,jk) = xkrgb(2,irgb) * e3t(ji,jj,jk,Kmm) ekr(ji,jj,jk) = xkrgb(3,irgb) * e3t(ji,jj,jk,Kmm) END_3D ! Heat flux across w-level (used in the dynamics) ! ----------------------------------------------- IF( ln_qsr_bio ) THEN ! zqsr_corr(:,:) = parsw(:,:) * qsr(:,:) ! ze0(:,:,1) = (1._wp - 3._wp * parsw(:,:)) * qsr(:,:) ! ( 1 - 3 * alpha ) * q ze1(:,:,1) = zqsr_corr(:,:) ze2(:,:,1) = zqsr_corr(:,:) ze3(:,:,1) = zqsr_corr(:,:) ! DO jk = 2, nksrp + 1 DO_2D(1, 1, 1, 1) ze0(ji,jj,jk) = ze0(ji,jj,jk-1) * EXP( -e3t(ji,jj,jk-1,Kmm) * (1. / rn_si0) ) ze1(ji,jj,jk) = ze1(ji,jj,jk-1) * EXP( -ekb (ji,jj,jk-1 ) ) ze2(ji,jj,jk) = ze2(ji,jj,jk-1) * EXP( -ekg (ji,jj,jk-1 ) ) ze3(ji,jj,jk) = ze3(ji,jj,jk-1) * EXP( -ekr (ji,jj,jk-1 ) ) END_2D END DO ! etot3(:,:,1) = qsr(:,:) * tmask(:,:,1) DO jk = 2, nksrp + 1 etot3(:,:,jk) = ( ze0(:,:,jk) + ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk) ) * tmask(:,:,jk) END DO ! ! ------------------------ ENDIF ! Photosynthetically Available Radiation (PAR) ! -------------------------------------------- zqsr_corr(:,:) = parsw(:,:) * qsr(:,:) / ( 1.-fr_i(:,:) + rtrn ) ! CALL trc_opt_par( kt, zqsr_corr, ze1, ze2, ze3 ) ! DO jk = 1, nksrp etot (:,:,jk) = ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk) ENDDO ! No Diurnal cycle PAR IF( l_trcdm2dc ) THEN zqsr_corr(:,:) = parsw(:,:) * qsr_mean(:,:) / ( 1.-fr_i(:,:) + rtrn ) ! CALL trc_opt_par( kt, zqsr_corr, ze1, ze2, ze3 ) DO jk = 1, nksrp etot_ndcy(:,:,jk) = ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk) END DO ELSE etot_ndcy(:,:,:) = etot(:,:,:) ENDIF ! Weighted broadband attenuation coefficient ! ------------------------------------------ DO_3D( 1, 1, 1, 1, 1, jpkm1 ) ztmp = ze1(ji,jj,jk)* ekb(ji,jj,jk) + ze2(ji,jj,jk) * ekg(ji,jj,jk) + ze3(ji,jj,jk) * ekr(ji,jj,jk) xeps(ji,jj,jk) = ztmp / e3t(ji,jj,jk,Kmm) / (etot(ji,jj,jk) + rtrn) END_3D !xeps = (ze1(:,:,:) * ekb(:,:,:) + ze2(:,:,:) * ekg(:,:,:) + ze3(:,:,:) * ekr(:,:,:)) / e3t(:,:,:,Kmm) / (etot(:,:,:) + rtrn) ! Light at the euphotic depth ! --------------------------- zqsr100 = 0.01 * 3. * zqsr_corr(:,:) ! Euphotic depth and level ! ------------------------ neln (:,:) = 1 heup (:,:) = gdepw(:,:,2,Kmm) heup_01(:,:) = gdepw(:,:,2,Kmm) ! DO_3D( 1, 1, 1, 1, 2, nksrp ) IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= zqsr100(ji,jj) ) THEN ! Euphotic level (1st T-level strictly below Euphotic layer) ! NOTE: ensure compatibility with nmld_trc definition in trdmxl_trc neln(ji,jj) = jk+1 ! ! Euphotic layer depth heup(ji,jj) = gdepw(ji,jj,jk+1,Kmm) ENDIF ! Euphotic layer depth (light level definition) IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= 0.50 ) THEN heup_01(ji,jj) = gdepw(ji,jj,jk+1,Kmm) ENDIF END_3D ! heup (:,:) = MIN( 300., heup (:,:) ) heup_01(:,:) = MIN( 300., heup_01(:,:) ) ! IF( lk_iomput ) THEN CALL iom_put( "xeps" , xeps(:,:,:) * tmask(:,:,:) ) CALL iom_put( "Heup" , heup(:,: ) * tmask(:,:,1) ) ENDIF ! IF( ln_timing ) CALL timing_stop('trc_opt') ! END SUBROUTINE trc_opt SUBROUTINE trc_opt_par( kt, zqsr, pe1, pe2, pe3) !!---------------------------------------------------------------------- !! *** routine trc_opt_par *** !! !! ** purpose : compute PAR of each wavelength (Red-Green-Blue) !! from given surface shortwave radiation !! !!---------------------------------------------------------------------- INTEGER , INTENT(in) :: kt ! ocean time-step REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: zqsr ! real shortwave REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(out) :: pe1 , pe2 , pe3 ! PAR (R-G-B) ! INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp), DIMENSION(jpi,jpj) :: we1, we2, we3 ! PAR (R-G-B) at w-level !!---------------------------------------------------------------------- pe1(:,:,:) = 0. ; pe2(:,:,:) = 0. ; pe3(:,:,:) = 0. ! IF ( TRIM(light_loc) == 'center' ) THEN ! cell-center (t-pivot) pe1(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekb(:,:,1) ) pe2(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekg(:,:,1) ) pe3(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekr(:,:,1) ) ! DO_3D( 1, 1, 1, 1, 2, nksrp ) pe1(ji,jj,jk) = pe1(ji,jj,jk-1) * EXP( -0.5 * ( ekb(ji,jj,jk-1) + ekb(ji,jj,jk) ) ) pe2(ji,jj,jk) = pe2(ji,jj,jk-1) * EXP( -0.5 * ( ekg(ji,jj,jk-1) + ekg(ji,jj,jk) ) ) pe3(ji,jj,jk) = pe3(ji,jj,jk-1) * EXP( -0.5 * ( ekr(ji,jj,jk-1) + ekr(ji,jj,jk) ) ) END_3D ! ELSE IF ( TRIM(light_loc) == 'integral' ) THEN ! integrate over cell thickness we1(:,:) = zqsr(:,:) we2(:,:) = zqsr(:,:) we3(:,:) = zqsr(:,:) ! DO_3D( 1, 1, 1, 1, 1, nksrp ) ! integrate PAR over current t-level pe1(ji,jj,jk) = we1(ji,jj) / (ekb(ji,jj,jk) + rtrn) * (1. - EXP( -ekb(ji,jj,jk) )) pe2(ji,jj,jk) = we2(ji,jj) / (ekg(ji,jj,jk) + rtrn) * (1. - EXP( -ekg(ji,jj,jk) )) pe3(ji,jj,jk) = we3(ji,jj) / (ekr(ji,jj,jk) + rtrn) * (1. - EXP( -ekr(ji,jj,jk) )) ! PAR at next w-level we1(ji,jj) = we1(ji,jj) * EXP( -ekb(ji,jj,jk) ) we2(ji,jj) = we2(ji,jj) * EXP( -ekg(ji,jj,jk) ) we3(ji,jj) = we3(ji,jj) * EXP( -ekr(ji,jj,jk) ) END_3D ! ENDIF ! ! END SUBROUTINE trc_opt_par SUBROUTINE trc_opt_sbc( kt ) !!---------------------------------------------------------------------- !! *** routine trc_opt_sbc *** !! !! ** purpose : read and interpolate the variable PAR fraction !! of shortwave radiation !! !! ** method : read the files and interpolate the appropriate variables !! !! ** input : external netcdf files !! !!---------------------------------------------------------------------- INTEGER, INTENT(in) :: kt ! ocean time step ! INTEGER :: ji,jj REAL(wp) :: zcoef !!--------------------------------------------------------------------- ! IF( ln_timing ) CALL timing_start('trc_opt_sbc') ! ! Compute par_varsw at nit000 or only if there is more than 1 time record in par coefficient file IF( ln_varpar ) THEN IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_par > 1 ) ) THEN CALL fld_read( kt, 1, sf_par ) par_varsw(:,:) = ( sf_par(1)%fnow(:,:,1) ) / 3.0 ENDIF ENDIF ! IF( ln_timing ) CALL timing_stop('trc_opt_sbc') ! END SUBROUTINE trc_opt_sbc SUBROUTINE trc_opt_ini !!---------------------------------------------------------------------- !! *** ROUTINE trc_opt_ini *** !! !! ** Purpose : Initialization of tabulated attenuation coefficients !! and percentage of PAR in Shortwave !! !! ** Input : external ascii and netcdf files !!---------------------------------------------------------------------- INTEGER :: numpar, ierr, ios ! Local integer ! CHARACTER(len=100) :: cn_dir ! Root directory for location of ssr files TYPE(FLD_N) :: sn_par ! informations about the fields to be read ! NAMELIST/namtrc_opt/cn_dir, sn_par, ln_varpar, parlux, light_loc !!---------------------------------------------------------------------- IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) 'trc_opt_ini : Initialize light module' WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' ENDIF READ ( numnat_ref, namtrc_opt, IOSTAT = ios, ERR = 901) 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_opt in reference namelist' ) READ ( numnat_cfg, namtrc_opt, IOSTAT = ios, ERR = 902 ) 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtrc_opt in configuration namelist' ) IF(lwm) WRITE ( numont, namtrc_opt ) IF(lwp) THEN WRITE(numout,*) ' Namelist : namtrc_opt ' WRITE(numout,*) ' PAR as a variable fraction of SW ln_varpar = ', ln_varpar WRITE(numout,*) ' Fraction of shortwave as PAR parlux = ', parlux WRITE(numout,*) ' Light location in the water cell light_loc = ', light_loc ENDIF ! ! Variable PAR at the surface of the ocean ! ---------------------------------------- IF( ln_varpar ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) ' ==>>> initialize variable par fraction (ln_varpar=T)' ! ALLOCATE( par_varsw(jpi,jpj) ) ! ALLOCATE( sf_par(1), STAT=ierr ) !* allocate and fill sf_par (forcing structure) with sn_par IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'trc_opt_ini: unable to allocate sf_par structure' ) ! CALL fld_fill( sf_par, (/ sn_par /), cn_dir, 'trc_opt_ini', 'Initialize prescribed PAR forcing ', 'namtrc_opt' ) ALLOCATE( sf_par(1)%fnow(jpi,jpj,1) ) IF( sn_par%ln_tint ) ALLOCATE( sf_par(1)%fdta(jpi,jpj,1,2) ) CALL iom_open ( TRIM( sn_par%clname ) , numpar ) ntimes_par = iom_getszuld( numpar ) ! get number of record in file ENDIF ! CALL trc_oce_rgb( xkrgb ) ! tabulated attenuation coefficients nksrp = trc_oce_ext_lev( r_si2, 0.33e2 ) ! max level of light extinction (Blue Chl=0.01) ! IF(lwp) WRITE(numout,*) ' level of light extinction = ', nksrp, ' ref depth = ', gdepw_1d(nksrp+1), ' m' ! ekr (:,:,:) = 0._wp ekb (:,:,:) = 0._wp ekg (:,:,:) = 0._wp etot (:,:,:) = 0._wp etot_ndcy(:,:,:) = 0._wp IF( ln_qsr_bio ) etot3 (:,:,:) = 0._wp ! END SUBROUTINE trc_opt_ini INTEGER FUNCTION trc_opt_alloc() !!---------------------------------------------------------------------- !! *** ROUTINE trc_opt_alloc *** !!---------------------------------------------------------------------- ! ALLOCATE( ekb(jpi,jpj,jpk), ekr(jpi,jpj,jpk), & ekg(jpi,jpj,jpk), xeps(jpi,jpj,jpk), STAT= trc_opt_alloc ) ! IF( trc_opt_alloc /= 0 ) CALL ctl_stop( 'STOP', 'trc_opt_alloc : failed to allocate arrays.' ) ! END FUNCTION trc_opt_alloc !!====================================================================== END MODULE trcopt