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.
trcopt_medusa.F90 in branches/UKMO/dev_r5518_GO6_package_FOAMv14_readchl/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

source: branches/UKMO/dev_r5518_GO6_package_FOAMv14_readchl/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcopt_medusa.F90 @ 14586

Last change on this file since 14586 was 14586, checked in by dford, 3 years ago

Allow just physics to use input chlorophyll.
File size: 10.8 KB
Line 
1MODULE trcopt_medusa
2   !!======================================================================
3   !!                         ***  MODULE trcopt_medusa  ***
4   !! TOP :   MEDUSA Compute the light availability in the water column
5   !!======================================================================
6   !! History :    -   !  1995-05  (M. Levy) Original code
7   !!              -   !  1999-09  (J.-M. Andre, M. Levy)
8   !!              -   !  1999-11  (C. Menkes, M.-A. Foujols) itabe initial
9   !!              -   !  2000-02  (M.A. Foujols) change x**y par exp(y*log(x))
10   !!             2.0  !  2007-12  (C. Deltel, G. Madec)  F90
11   !!              -   !  2008-08  (K. Popova) adaptation for MEDUSA
12   !!              -   !  2008-11  (A. Yool) continuing adaptation for MEDUSA
13   !!              -   !  2010-03  (A. Yool) updated for branch inclusion
14   !!----------------------------------------------------------------------
15#if defined key_medusa
16   !!----------------------------------------------------------------------
17   !!   'key_medusa'                                      MEDUSA bio-model
18   !!----------------------------------------------------------------------
19   !!   trc_opt_medusa        :   Compute the light availability in the water column
20   !!----------------------------------------------------------------------
21   USE oce_trc         !
22   USE trc
23   USE prtctl_trc      ! Print control for debbuging
24   USE sms_medusa
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   trc_opt_medusa   ! called in trcprg.F90
30
31   REAL(wp), DIMENSION(3,61) :: okrgb   !: tabulated attenuation coefficients for RGB absorption
32
33   !!* Substitution
34#  include "domzgr_substitute.h90"
35   !!----------------------------------------------------------------------
36   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007)
37   !! $Id$
38   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
39   !!----------------------------------------------------------------------
40
41CONTAINS
42
43   SUBROUTINE trc_opt_medusa( kt )
44      !!---------------------------------------------------------------------
45      !!                     ***  ROUTINE trc_opt_medusa  ***
46      !!
47      !! ** Purpose :   computes the light propagation in the water column
48      !!              and the euphotic layer depth
49      !!
50      !! ** Method  :   local par is computed in w layers using light propagation
51      !!              mean par in t layers are computed by integration
52      !!---------------------------------------------------------------------
53      INTEGER, INTENT( in ) ::   kt   ! index of the time stepping
54      INTEGER  ::   ji, jj, jk, irgb
55      REAL(wp) ::   zpig                                    ! total pigment
56      REAL(wp) ::   zkr                                     ! total absorption coefficient in red
57      REAL(wp) ::   zkg                                     ! total absorption coefficient in green
58      REAL(wp) ::   totchl                                  ! total Chl concentreation
59      REAL(wp), DIMENSION(jpi,jpj)     ::   zpar100         ! irradiance at euphotic layer depth
60      REAL(wp), DIMENSION(jpi,jpj)     ::   zpar0m          ! irradiance just below the surface
61      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zparr, zparg    ! red and green compound of par
62      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zekb, zekg, zekr
63      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze1, ze2, ze3
64
65      CHARACTER (len=25) :: charout
66      !!---------------------------------------------------------------------
67
68      !! AXY (20/11/14): alter this to report on first MEDUSA call
69      !! IF( kt == nit000 ) THEN
70      IF( kt == nittrc000 ) THEN
71         IF(lwp) WRITE(numout,*)
72         IF(lwp) WRITE(numout,*) ' trc_opt_medusa: MEDUSA optic-model'
73         IF(lwp) WRITE(numout,*) ' ~~~~~~~'
74    IF(lwp) WRITE(numout,*) ' kt =',kt
75      ENDIF
76
77      ! determination of surface irradiance
78      ! -----------------------------------
79      ! AXY (23/07/15): the inclusion of empirical DMS calculations requires
80      !                 daily averages of a series of properties that are
81      !                 used as inputs; these include surface irradiance;
82      !                 here, this is taken advantage of to allow MEDUSA to
83      !                 base its submarine light field on daily average
84      !                 rather than "instantaneous" irradiance; largely
85      !                 because MEDUSA was originally formulated to work
86      !                 with diel average irradiance rather than a diel
87      !                 cycle; using key_avgqsr_medusa activates this
88      !                 functionality, while its absence gives default
89      !                 MEDUSA (which is whatever is supplied by NEMO)
90# if defined key_avgqsr_medusa
91      ! surface irradiance input is rolling average irradiance
92      zpar0m (:,:)   = zn_dms_qsr(:,:) * 0.43
93# else     
94      ! surface irradiance input is   instantaneous irradiance
95      zpar0m (:,:)   =        qsr(:,:) * 0.43
96# endif
97      ! AXY (22/08/14): when zpar0m = 0, zpar100 is also zero and calculating
98      !                 euphotic depth is not possible (cf. the Arctic Octopus);
99      !                 a "solution" to this is to set zpar0m to some minimal
100      !                 value such that zpar100 also has a non-zero value and
101      !                 euphotic depth can be calculated properly; note that,
102      !                 in older, non-diurnal versions of NEMO, this was much
103      !                 less of a problem; note also that, if pushed, I will
104      !                 claim that my minimal value of zpar0m refers to light
105      !                 from stars
106      DO jj = 1, jpj
107         DO ji = 1, jpi
108            IF( zpar0m(ji,jj) <= 0.0 ) zpar0m(ji,jj) = 0.001  ! = 1 mW/m2
109         ENDDO
110      ENDDO
111      zpar100(:,:)   = zpar0m(:,:) * 0.01
112      xpar   (:,:,1) = zpar0m(:,:)
113      zparr  (:,:,1) = 0.5 * zpar0m(:,:)
114      zparg  (:,:,1) = 0.5 * zpar0m(:,:)
115
116      ! determination of xpar
117      ! ---------------------
118
119      IF ( ln_rgb ) THEN
120         ! Mean PAR in T levels using RGB scheme
121         ! Should be same as in traqsr, but T rather than W levels
122         IF( kt == nittrc000 ) THEN
123            CALL trc_oce_rgb( okrgb )
124         ENDIF
125         DO jk = 1, jpkm1
126            DO jj = 1, jpj
127               DO ji = 1, jpi
128                  IF ( ln_chlbio ) THEN
129                     totchl = chl_for_qsr(ji,jj,jk)
130                  ELSE
131                     totchl = trn(ji,jj,jk,jpchn) + trn(ji,jj,jk,jpchd)
132                  ENDIF
133                  totchl = MIN( 10. , MAX( 0.05, totchl ) )
134                  irgb = NINT( 41 + 20.* LOG10( totchl ) + 1.e-15 )
135                  !                                                         
136                  zekb(ji,jj,jk) = okrgb(1,irgb) * fse3t(ji,jj,jk)
137                  zekg(ji,jj,jk) = okrgb(2,irgb) * fse3t(ji,jj,jk)
138                  zekr(ji,jj,jk) = okrgb(3,irgb) * fse3t(ji,jj,jk)
139               END DO
140            END DO
141         END DO
142         ze1(:,:,1) = zpar0m(:,:) * EXP( -0.5 * zekb(:,:,1) ) / 3.0
143         ze2(:,:,1) = zpar0m(:,:) * EXP( -0.5 * zekg(:,:,1) ) / 3.0
144         ze3(:,:,1) = zpar0m(:,:) * EXP( -0.5 * zekr(:,:,1) ) / 3.0
145         !
146         DO jk = 2, jpk
147            DO jj = 1, jpj
148               DO ji = 1, jpi
149                  ze1(ji,jj,jk) = ze1(ji,jj,jk-1) * EXP( -0.5 * ( zekb(ji,jj,jk-1) + zekb(ji,jj,jk) ) )
150                  ze2(ji,jj,jk) = ze2(ji,jj,jk-1) * EXP( -0.5 * ( zekg(ji,jj,jk-1) + zekg(ji,jj,jk) ) )
151                  ze3(ji,jj,jk) = ze3(ji,jj,jk-1) * EXP( -0.5 * ( zekr(ji,jj,jk-1) + zekr(ji,jj,jk) ) )
152                  xpar(ji,jj,jk) = MAX( ze1(ji,jj,jk) + ze2(ji,jj,jk) + ze3(ji,jj,jk), 1.e-15 )
153               END DO
154            END DO
155         END DO
156      ELSE
157         DO jk = 2, jpk                     ! determination of local par in w levels
158            DO jj = 1, jpj
159               DO ji = 1, jpi
160                  totchl =trn(ji,jj,jk-1,jpchn)+trn(ji,jj,jk-1,jpchd)
161                  zpig = MAX( TINY(0.), totchl/rpig) 
162                  zkr  = xkr0 + xkrp * EXP( xlr * LOG( zpig ) )
163                  zkg  = xkg0 + xkgp * EXP( xlg * LOG( zpig ) )
164                  zparr(ji,jj,jk) = zparr(ji,jj,jk-1) * EXP( -zkr * fse3t(ji,jj,jk-1) )
165                  zparg(ji,jj,jk) = zparg(ji,jj,jk-1) * EXP( -zkg * fse3t(ji,jj,jk-1) )
166               END DO
167           END DO
168         END DO
169
170         DO jk = 1, jpkm1                   ! mean par in t levels
171            DO jj = 1, jpj
172               DO ji = 1, jpi
173                  totchl =trn(ji,jj,jk  ,jpchn)+trn(ji,jj,jk  ,jpchd)
174                  zpig = MAX( TINY(0.), totchl/rpig) 
175                  zkr  = xkr0 + xkrp * EXP( xlr * LOG( zpig ) )
176                  zkg  = xkg0 + xkgp * EXP( xlg * LOG( zpig ) )
177                  zparr(ji,jj,jk)    = zparr(ji,jj,jk) / zkr / fse3t(ji,jj,jk) * ( 1 - EXP( -zkr*fse3t(ji,jj,jk) ) )
178                  zparg(ji,jj,jk)    = zparg(ji,jj,jk) / zkg / fse3t(ji,jj,jk) * ( 1 - EXP( -zkg*fse3t(ji,jj,jk) ) )
179                  xpar (ji,jj,jk) = MAX( zparr(ji,jj,jk) + zparg(ji,jj,jk), 1.e-15 )
180               END DO
181            END DO
182         END DO
183      ENDIF
184
185      ! 3. Determination of euphotic layer depth
186      ! ----------------------------------------
187
188      ! Euphotic layer bottom level
189      neln(:,:) = 1                                           ! initialisation of EL level
190      DO jk = 1, jpkm1
191         DO jj = 1, jpj
192           DO ji = 1, jpi
193              IF( xpar(ji,jj,jk) >= zpar100(ji,jj) )   neln(ji,jj) = jk+1 ! 1rst T-level strictly below EL bottom
194              !                                                  ! nb. this is to ensure compatibility with
195              !                                                  ! nmld_trc definition in trd_mld_trc_zint
196           END DO
197         END DO
198      ENDDO
199
200      ! Euphotic layer depth
201      !! Jpalm -- 06-03-2017 -- add init xze, to avoid halo problems within the
202      !!                        writing process
203      xze(:,:) = 0.0
204      DO jj = 1, jpj
205         DO ji = 1, jpi
206            xze(ji,jj) = fsdepw( ji, jj, neln(ji,jj) )            ! exact EL depth
207         END DO
208      ENDDO 
209
210      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
211         WRITE(charout, FMT="('opt')")
212         CALL prt_ctl_trc_info(charout)
213         CALL prt_ctl_trc( tab4d=trn, mask=tmask, clinfo=ctrcnm )
214      ENDIF
215
216   END SUBROUTINE trc_opt_medusa
217
218#else
219   !!======================================================================
220   !!  Dummy module :                                   No MEDUSA bio-model
221   !!======================================================================
222CONTAINS
223   SUBROUTINE trc_opt_medusa( kt )                   ! Empty routine
224      INTEGER, INTENT( in ) ::   kt
225      WRITE(*,*) 'trc_opt_medusa: You should not have seen this print! error?', kt
226   END SUBROUTINE trc_opt_medusa
227#endif 
228
229   !!======================================================================
230END MODULE  trcopt_medusa
Note: See TracBrowser for help on using the repository browser.