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_under_ice_relax_dr_hook/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcopt_medusa.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 5 years ago

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