Changeset 14415 for NEMO/branches/2020/dev_r12985_TOP-04_IMMERSE_BGC_interface/src/TOP/trcopt.F90

Timestamp:

2021-02-06T14:04:27+01:00 (4 years ago)

Author:

lovato

Message:

trcopt revised computation of light at different location (ticket #2489)

File:

• NEMO/branches/2020/dev_r12985_TOP-04_IMMERSE_BGC_interface/src/TOP/trcopt.F90 (modified) (5 diffs)

NEMO/branches/2020/dev_r12985_TOP-04_IMMERSE_BGC_interface/src/TOP/trcopt.F90

@@ -24 +24 @@
    LOGICAL  ::   ln_varpar   ! boolean for variable PAR fraction
    REAL(wp) ::   parlux      ! Fraction of shortwave as PAR
-   CHARACTER (len=25) :: light_loc
+   CHARACTER (len=25) :: light_loc ! Light location in the water cell ('center', 'integral')
 
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_par        ! structure of input par
@@ -133 +133 @@
       ENDDO
 
-      ! weighted broadband diffuse attenuation coefficient
-      WHERE(etot .ne. 0.) &
-         xeps = ( ze1 * ekr + ze2 * ekb + ze3 * ekg ) / e3t(:,:,:,Kmm) / etot
-      ! Compute PAR at cell center (T-level) or integrate over cell depth
-      IF ( TRIM(light_loc) == 'center' ) THEN
-          IF(lwp) WRITE(numout,*) 'trcopt : center'
-          etot = etot * EXP( -xeps * 0.5 * e3t(:,:,:,Kmm))
-      ELSE IF ( TRIM(light_loc) == 'integral' ) THEN
-          IF(lwp) WRITE(numout,*) 'trcopt : integral'
-          WHERE(etot == 0.) &
-          etot = etot / xeps * (1. - EXP(-xeps*e3t(:,:,:,Kmm)))
-      ENDIF
-
-      !
       ! No Diurnal cycle PAR
       IF( l_trcdm2dc ) THEN
@@ -158 +144 @@
          etot_ndcy(:,:,:) = etot(:,:,:)
       ENDIF
+
+      !     Weighted broadband attenuation coefficient
+      !     ------------------------------------------
+      xeps = (ze1(:,:,:) * ekb(:,:,:) + ze2(:,:,:) * ekg(:,:,:) + ze3(:,:,:) * ekr(:,:,:)) / e3t(:,:,:,Kmm) / (etot(:,:,:) + rtrn)
 
       !     Light at the euphotic depth
@@ -202 +192 @@
       !!
       !! ** purpose :   compute PAR of each wavelength (Red-Green-Blue)
-      !!                for a given shortwave radiation at w-level
+      !!                from given surface shortwave radiation
       !!
       !!----------------------------------------------------------------------
@@ -209 +199 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(out)     ::   pe1 , pe2 , pe3   ! PAR (R-G-B)
       !
-      INTEGER    ::   ji, jj, jk                   ! dummy loop indices
+      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.
       !
-      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
+      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 
       !
       !