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 @ 230

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

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