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.
traqsr.F90 in trunk/NEMO/OPA_SRC/TRA – NEMO

source: trunk/NEMO/OPA_SRC/TRA/traqsr.F90 @ 247

Last change on this file since 247 was 247, checked in by opalod, 19 years ago

CL : Add CVS Header and CeCILL licence information
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 15.7 KB
Line 
1MODULE traqsr
2   !!======================================================================
3   !!                       ***  MODULE  traqsr  ***
4   !! Ocean physics: solar radiation penetration in the top ocean levels
5   !!======================================================================
6
7   !!----------------------------------------------------------------------
8   !!   tra_qsr      : trend due to the solar radiation penetration
9   !!   tra_qsr_init : solar radiation penetration initialization
10   !!----------------------------------------------------------------------
11   !! * Modules used
12   USE oce             ! ocean dynamics and active tracers
13   USE dom_oce         ! ocean space and time domain
14   USE trdmod          ! ocean active tracers trends
15   USE trdmod_oce      ! ocean variables trends
16   USE in_out_manager  ! I/O manager
17
18   USE trc_oce         ! share SMS/Ocean variables
19
20   USE ocesbc          ! thermohaline fluxes
21   USE phycst          ! physical constants
22
23   IMPLICIT NONE
24   PRIVATE
25
26   !! * Routine accessibility
27   PUBLIC tra_qsr      ! routine called by step.F90 (ln_traqsr=T)
28   PUBLIC tra_qsr_init ! routine called by opa.F90
29
30   !! * Shared module variables
31   LOGICAL, PUBLIC ::   ln_traqsr = .TRUE.   !: qsr flag (Default=T)
32
33   !! * Module variables
34   REAL(wp), PUBLIC ::  & !!! * penetrative solar radiation namelist *
35      rabs = 0.58_wp,   &  ! fraction associated with xsi1
36      xsi1 = 0.35_wp,   &  ! first depth of extinction
37      xsi2 = 23.0_wp       ! second depth of extinction
38      !                    ! (default values: water type Ib)
39   LOGICAL ::           & 
40      ln_qsr_sms = .false. ! flag to use or not the biological
41      !                    ! fluxes for light
42   
43   INTEGER ::    &
44      nksr                 ! number of levels
45   REAL(wp), DIMENSION(jpk) ::   &
46      gdsr                 ! profile of the solar flux penetration
47
48   !! * Substitutions
49#  include "domzgr_substitute.h90"
50#  include "vectopt_loop_substitute.h90"
51   !!----------------------------------------------------------------------
52   !!   OPA 9.0 , LOCEAN-IPSL (2005)
53   !! $Header$
54   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
55   !!----------------------------------------------------------------------
56
57CONTAINS
58
59   SUBROUTINE tra_qsr( kt )
60      !!----------------------------------------------------------------------
61      !!                  ***  ROUTINE tra_qsr  ***
62      !!
63      !! ** Purpose :   Compute the temperature trend due to the solar radiation
64      !!      penetration and add it to the general temperature trend.
65      !!
66      !! ** Method  : The profile of the solar radiation within the ocean is
67      !!      defined through two penetration length scale (xsr1,xsr2) and a
68      !!      ratio (rabs) as :
69      !!         I(k) = Qsr*( rabs*EXP(z(k)/xsr1) + (1.-rabs)*EXP(z(k)/xsr2) )
70      !!         The temperature trend associated with the solar radiation
71      !!      penetration is given by :
72      !!            zta = 1/e3t dk[ I ] / (rau0*Cp)
73      !!         At the bottom, boudary condition for the radiation is no flux :
74      !!      all heat which has not been absorbed in the above levels is put
75      !!      in the last ocean level.
76      !!         In z-coordinate case, the computation is only done down to the
77      !!      level where I(k) < 1.e-15 W/m2. In addition, the coefficients
78      !!      used for the computation are calculated one for once as they
79      !!      depends on k only.
80      !!
81      !! ** Action  : - update ta with the penetrative solar radiation trend
82      !!              - save the trend in ttrd ('key_trdtra')
83      !!
84      !! History :
85      !!   6.0  !  90-10  (B. Blanke)  Original code
86      !!   7.0  !  91-11  (G. Madec)
87      !!        !  96-01  (G. Madec)  s-coordinates
88      !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
89      !!   9.0  !  04-08  (C. Talandier) New trends organization
90      !!----------------------------------------------------------------------
91      !! * Modules used     
92      USE oce, ONLY :    ztdta => ua,   & ! use ua as 3D workspace   
93                         ztdsa => va      ! use va as 3D workspace   
94
95      !! * Arguments
96      INTEGER, INTENT( in ) ::   kt       ! ocean time-step
97
98      !! * Local declarations
99      INTEGER ::    ji, jj, jk            ! dummy loop indexes
100      REAL(wp) ::   zc0, zta              ! temporary scalars
101      REAL(wp) ::   zc1 , zc2 ,        &  ! temporary scalars
102                    zdp1, zdp2            !
103      !!----------------------------------------------------------------------
104
105      IF( kt == nit000 ) THEN
106         IF ( lwp )  WRITE(numout,*)
107         IF ( lwp )  WRITE(numout,*) 'tra_qsr : penetration of the surface solar radiation'
108         IF ( lwp )  WRITE(numout,*) '~~~~~~~'
109      ENDIF
110
111      ! Save ta and sa trends
112      IF( l_trdtra )   THEN
113         ztdta(:,:,:) = ta(:,:,:) 
114         ztdsa(:,:,:) = 0.e0
115      ENDIF
116
117      IF( lk_qsr_sms .AND. ln_qsr_sms ) THEN    !  Biological fluxes  !
118         !                                      ! =================== !
119         !
120         !                                                ! ===============
121         DO jk = 1, jpkm1                                 ! Horizontal slab
122            !                                             ! ===============
123            DO jj = 2, jpjm1
124               DO ji = fs_2, fs_jpim1   ! vector opt.
125                 
126                  zc0 = ro0cpr  / fse3t(ji,jj,jk)      ! compute the qsr trend
127                  zta = zc0 * ( etot3(ji,jj,jk  ) * tmask(ji,jj,jk)     &
128                     &        - etot3(ji,jj,jk+1) * tmask(ji,jj,jk+1) )
129                 
130                  ta(ji,jj,jk) = ta(ji,jj,jk) + zta       ! add qsr trend to the temperature trend
131                 
132               END DO
133            END DO
134            !                                             ! ===============
135         END DO                                           !   End of slab
136         !                                                ! ===============
137         ! save the trends for diagnostic
138         ! qsr tracers trends
139         IF( l_trdtra )   THEN
140            ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:)
141            CALL trd_mod(ztdta, ztdsa, jpttdqsr, 'TRA', kt)
142         ENDIF
143
144      ELSE
145         !                                                ! =================== !
146         IF( lk_sco ) THEN                                !    s-coordinate     !
147         !                                                ! =================== !
148         !
149         !                                                   ! ===============
150            DO jk = 1, jpkm1                                 ! Horizontal slab
151            !                                                ! ===============
152               DO jj = 2, jpjm1
153                  DO ji = fs_2, fs_jpim1   ! vector opt.
154
155                     zdp1 = -fsdepw(ji,jj,jk  )              ! compute the qsr trend
156                     zdp2 = -fsdepw(ji,jj,jk+1)
157                     zc0  = qsr(ji,jj) * ro0cpr / fse3t(ji,jj,jk)
158                     zc1  =   (  rabs * EXP(zdp1/xsi1) + (1.-rabs) * EXP(zdp1/xsi2)  )
159                     zc2  = - (  rabs * EXP(zdp2/xsi1) + (1.-rabs) * EXP(zdp2/xsi2)  )
160                     zta  = zc0 * (  zc1 * tmask(ji,jj,jk) + zc2 * tmask(ji,jj,jk+1)  )
161                     
162                     ta(ji,jj,jk) = ta(ji,jj,jk) + zta       ! add qsr trend to the temperature trend
163                     
164                  END DO
165               END DO
166               !                                                ! ===============
167            END DO                                              !   End of slab
168            !                                                   ! ===============
169            ! save the trends for diagnostic
170            ! qsr tracers trends
171            IF( l_trdtra )   THEN
172               ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:)
173               CALL trd_mod(ztdta, ztdsa, jpttdqsr, 'TRA', kt)
174            ENDIF
175            !
176         ENDIF
177         !                                                ! =================== !
178         IF( lk_zps ) THEN                                !    partial steps    !
179            !                                             ! =================== !
180            !
181            !                                                ! ===============
182            DO jk = 1, nksr                                  ! Horizontal slab
183               !                                             ! ===============
184               DO jj = 2, jpjm1
185                  DO ji = fs_2, fs_jpim1   ! vector opt.
186                     
187                     zc0 = qsr(ji,jj) / fse3t(ji,jj,jk)      ! compute the qsr trend
188                     zta = zc0 * ( gdsr(jk) * tmask(ji,jj,jk) - gdsr(jk+1) * tmask(ji,jj,jk+1) )
189                     
190                     ta(ji,jj,jk) = ta(ji,jj,jk) + zta       ! add qsr trend to the temperature trend
191                     
192                  END DO
193               END DO
194               !                                             ! ===============
195            END DO                                           !   End of slab
196            !                                                ! ===============
197            ! save the trends for diagnostic
198            ! qsr tracers trends
199            IF( l_trdtra )   THEN
200               ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:)
201               CALL trd_mod(ztdta, ztdsa, jpttdqsr, 'TRA', kt)
202            ENDIF
203            !
204         ENDIF
205         !                                                ! =================== !
206         IF( lk_zco ) THEN                                !     z-coordinate    !
207            !                                             ! =================== !
208            !
209            !                                                ! ===============
210            DO jk = 1, nksr                                  ! Horizontal slab
211               !                                             ! ===============
212               zc0 = 1. / fse3t(1,1,jk)
213               DO jj = 2, jpjm1
214                  DO ji = fs_2, fs_jpim1   ! vector opt.
215                     !                                       ! compute qsr forcing trend
216                     zta = qsr(ji,jj) * zc0 * ( gdsr(jk)*tmask(ji,jj,jk) - gdsr(jk+1)*tmask(ji,jj,jk+1) )
217                     
218                     ta(ji,jj,jk) = ta(ji,jj,jk) + zta       ! add qsr trend to the temperature trend
219                     
220                  END DO
221               END DO
222               !                                             ! ===============
223            END DO                                           !   End of slab
224            !                                                ! ===============
225            ! save the trends for diagnostic
226            ! qsr tracers trends
227            IF( l_trdtra )   THEN
228               ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:)
229               CALL trd_mod(ztdta, ztdsa, jpttdqsr, 'TRA', kt)
230            ENDIF
231            !
232         ENDIF
233         !
234      ENDIF
235
236
237      IF(l_ctl) THEN         ! print mean trends (used for debugging)
238         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
239         WRITE(numout,*) ' qsr  - Ta: ', zta-t_ctl
240         t_ctl = zta 
241      ENDIF
242
243
244      ! 2. Substract qsr from the total heat flux q
245      ! -------------------------------------------
246      DO jj = 2, jpjm1
247         DO ji = fs_2, fs_jpim1   ! vector opt.
248            q(ji,jj) = qt(ji,jj) - qsr(ji,jj)
249         END DO
250      END DO
251
252   END SUBROUTINE tra_qsr
253
254
255   SUBROUTINE tra_qsr_init
256      !!----------------------------------------------------------------------
257      !!                  ***  ROUTINE tra_qsr_init  ***
258      !!
259      !! ** Purpose :   Initialization for the penetrative solar radiation
260      !!
261      !! ** Method  :   The profile of solar radiation within the ocean is set
262      !!      from two length scale of penetration (xsr1,xsr2) and a ratio
263      !!      (rabs). These parameters are read in the namqsr namelist. The
264      !!      default values correspond to clear water (type I in Jerlov'
265      !!      (1968) classification.
266      !!         called by tra_qsr at the first timestep (nit000)
267      !!
268      !! ** Action  : - initialize xsr1, xsr2 and rabs
269      !!
270      !! Reference :
271      !!   Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
272      !!
273      !! History :
274      !!   8.5  !  02-06  (G. Madec) Original code
275      !!----------------------------------------------------------------------
276      !! * Local declarations
277      INTEGER ::    ji,jj,jk, &  ! dummy loop index
278                    indic        ! temporary integer
279      REAL(wp) ::   zdp1         ! temporary scalar
280
281      NAMELIST/namqsr/ ln_traqsr, rabs, xsi1, xsi2, ln_qsr_sms
282      !!----------------------------------------------------------------------
283
284      ! Read Namelist namqsr : ratio and length of penetration
285      ! --------------------
286      REWIND ( numnam )
287      READ   ( numnam, namqsr )
288
289      ! Parameter control and print
290      ! ---------------------------
291      IF( ln_traqsr  ) THEN
292        IF ( lwp ) THEN
293         WRITE(numout,*)
294         WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation'
295         WRITE(numout,*) '~~~~~~~~~~~~'
296         WRITE(numout,*) '    Namelist namqsr : set the parameter of penetration'
297         WRITE(numout,*) '        fraction associated with xsi     rabs        = ',rabs
298         WRITE(numout,*) '        first depth of extinction        xsi1        = ',xsi1
299         WRITE(numout,*) '        second depth of extinction       xsi2        = ',xsi2
300         IF( lk_qsr_sms ) THEN
301            WRITE(numout,*) '     Biological fluxes for light(Y/N) ln_qsr_sms  = ',ln_qsr_sms
302         ENDIF
303         WRITE(numout,*) ' '
304        END IF
305      ELSE
306        IF ( lwp ) THEN
307         WRITE(numout,*)
308         WRITE(numout,*) 'tra_qsr_init : NO solar flux penetration'
309         WRITE(numout,*) '~~~~~~~~~~~~'
310        END IF
311      ENDIF
312
313      IF( rabs > 1.e0 .OR. rabs < 0.e0 .OR. xsi1 < 0.e0 .OR. xsi2 < 0.e0 ) THEN
314         IF(lwp) WRITE(numout,cform_err)
315         IF(lwp) WRITE(numout,*) '             0<rabs<1, 0<xsi1, or 0<xsi2 not satisfied'
316         nstop = nstop + 1
317      ENDIF
318
319
320      ! Initialization
321      ! --------------
322      IF( .NOT. lk_sco ) THEN
323         ! z-coordinate with or without partial step : same before last ocean w-level everywhere
324         gdsr(:) = 0.e0
325         DO jk = 1, jpk
326            zdp1 = -fsdepw(1,1,jk)
327            gdsr(jk) = ro0cpr * (  rabs  * EXP( zdp1/xsi1 ) + (1.-rabs) * EXP( zdp1/xsi2 )  )
328            IF ( gdsr(jk) <= 1.e-10 ) EXIT
329         END DO
330         indic = 0
331         DO jk = 1, jpk
332            IF( gdsr(jk) <= 1.e-15 .AND. indic == 0 ) THEN
333               gdsr(jk) = 0.e0
334               nksr = jk
335               !!bug Edmee chg res   nksr = jk - 1
336               indic = 1
337            ENDIF
338         END DO
339         nksr = MIN( nksr, jpkm1 )
340         IF(lwp) THEN
341            WRITE(numout,*)
342            WRITE(numout,*) '        - z-coordinate, level max of computation =', nksr
343            WRITE(numout,*) '             profile of coef. of penetration:'
344            WRITE(numout,"('              ',7e11.2)") ( gdsr(jk), jk = 1, nksr )
345            WRITE(numout,*)
346         ENDIF
347         ! Initialisation of Biological fluxes for light here because
348         ! the optical biological model is call after the dynamical one
349         IF( lk_qsr_sms .AND. ln_qsr_sms ) THEN
350            DO jk = 1, jpkm1
351               DO jj = 1, jpj
352                  DO ji = 1, jpi
353                     etot3(ji,jj,jk) = qsr(ji,jj) * gdsr(jk) * tmask(ji,jj,jk) / ro0cpr
354                  END DO
355               END DO
356            END DO
357         ENDIF
358
359      ENDIF
360
361   END SUBROUTINE tra_qsr_init
362
363   !!======================================================================
364END MODULE traqsr
Note: See TracBrowser for help on using the repository browser.