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

Last change on this file since 72 was 32, checked in by opalod, 20 years ago

CT : UPDATE001 : First major NEMO update
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.4 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 trdtra_oce     ! ocean active tracer trends
15   USE in_out_manager  ! I/O manager
16
17   USE ocesbc          ! thermohaline fluxes
18   USE phycst          ! physical constants
19
20   IMPLICIT NONE
21   PRIVATE
22
23   !! * Routine accessibility
24   PUBLIC tra_qsr      ! routine called by step.F90 (ln_traqsr=T)
25   PUBLIC tra_qsr_init ! routine called by opa.F90
26
27   !! * Shared module variables
28   LOGICAL, PUBLIC ::   ln_traqsr = .TRUE.   !: qsr flag (Default=T)
29
30   !! * Module variables
31   REAL(wp) ::          & !!! * penetrative solar radiation namelist *
32      rabs = 0.58_wp,   &  ! fraction associated with xsi1
33      xsi1 = 0.35_wp,   &  ! first depth of extinction
34      xsi2 = 23.0_wp       ! second depth of extinction
35      !                    ! (default values: water type Ib)
36   INTEGER ::    &
37      nksr                 ! number of levels
38   REAL(wp), DIMENSION(jpk) ::   &
39      gdsr                 ! profile of the solar flux penetration
40
41   !! * Substitutions
42#  include "domzgr_substitute.h90"
43#  include "vectopt_loop_substitute.h90"
44   !!----------------------------------------------------------------------
45   !!   OPA 9.0 , LODYC-IPSL (2003)
46   !!----------------------------------------------------------------------
47
48CONTAINS
49
50   SUBROUTINE tra_qsr( kt )
51      !!----------------------------------------------------------------------
52      !!                  ***  ROUTINE tra_qsr  ***
53      !!
54      !! ** Purpose :   Compute the temperature trend due to the solar radiation
55      !!      penetration and add it to the general temperature trend.
56      !!
57      !! ** Method  : The profile of the solar radiation within the ocean is
58      !!      defined through two penetration length scale (xsr1,xsr2) and a
59      !!      ratio (rabs) as :
60      !!         I(k) = Qsr*( rabs*EXP(z(k)/xsr1) + (1.-rabs)*EXP(z(k)/xsr2) )
61      !!         The temperature trend associated with the solar radiation
62      !!      penetration is given by :
63      !!            zta = 1/e3t dk[ I ] / (rau0*Cp)
64      !!         At the bottom, boudary condition for the radiation is no flux :
65      !!      all heat which has not been absorbed in the above levels is put
66      !!      in the last ocean level.
67      !!         In z-coordinate case, the computation is only done down to the
68      !!      level where I(k) < 1.e-15 W/m2. In addition, the coefficients
69      !!      used for the computation are calculated one for once as they
70      !!      depends on k only.
71      !!
72      !! ** Action  : - update ta with the penetrative solar radiation trend
73      !!              - save the trend in ttrd ('key_trdtra')
74      !!
75      !! History :
76      !!   6.0  !  90-10  (B. Blanke)  Original code
77      !!   7.0  !  91-11  (G. Madec)
78      !!        !  96-01  (G. Madec)  s-coordinates
79      !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
80      !!----------------------------------------------------------------------
81      !! * Arguments
82      INTEGER, INTENT( in ) ::   kt    ! ocean time-step
83
84      !! * Local declarations
85      INTEGER ::    ji, jj, jk         ! dummy loop indexes
86      REAL(wp) ::   zc0, zta           ! temporary scalars
87      REAL(wp) ::   zc1 , zc2 ,     &  ! temporary scalars
88                    zdp1, zdp2         !
89      !!----------------------------------------------------------------------
90
91      IF( kt == nit000 ) THEN
92         IF ( lwp )  WRITE(numout,*)
93         IF ( lwp )  WRITE(numout,*) 'tra_qsr : penetration of the surface solar radiation'
94         IF ( lwp )  WRITE(numout,*) '~~~~~~~'
95      ENDIF
96
97      !                                                ! =================== !
98      IF( lk_sco ) THEN                                !     s-coordinate    !
99         !                                             ! =================== !
100         !
101         !                                                   ! ===============
102         DO jk = 1, jpkm1                                    ! Horizontal slab
103            !                                                ! ===============
104            DO jj = 2, jpjm1
105               DO ji = fs_2, fs_jpim1   ! vector opt.
106
107                  zdp1 = -fsdepw(ji,jj,jk  )              ! compute the qsr trend
108                  zdp2 = -fsdepw(ji,jj,jk+1)
109                  zc0  = qsr(ji,jj) * ro0cpr / fse3t(ji,jj,jk)
110                  zc1  =   (  rabs * EXP(zdp1/xsi1) + (1.-rabs) * EXP(zdp1/xsi2)  )
111                  zc2  = - (  rabs * EXP(zdp2/xsi1) + (1.-rabs) * EXP(zdp2/xsi2)  )
112                  zta  = zc0 * (  zc1 * tmask(ji,jj,jk) + zc2 * tmask(ji,jj,jk+1)  )
113
114                  ta(ji,jj,jk) = ta(ji,jj,jk) + zta       ! add qsr trend to the temperature trend
115
116#      if defined key_trdtra || defined key_trdmld
117                  ttrd(ji,jj,jk,7) = zta                  ! save the qsr trend
118#      endif
119               END DO
120            END DO
121            !                                                ! ===============
122         END DO                                              !   End of slab
123         !                                                   ! ===============
124      ENDIF
125      !                                                ! =================== !
126      IF( lk_zps ) THEN                                !    partial steps    !
127         !                                             ! =================== !
128         !
129         !                                                ! ===============
130         DO jk = 1, nksr                                  ! Horizontal slab
131            !                                             ! ===============
132            DO jj = 2, jpjm1
133               DO ji = fs_2, fs_jpim1   ! vector opt.
134   
135                  zc0 = qsr(ji,jj) / fse3t(ji,jj,jk)      ! compute the qsr trend
136                  zta = zc0 * ( gdsr(jk) * tmask(ji,jj,jk) - gdsr(jk+1) * tmask(ji,jj,jk+1) )
137   
138                  ta(ji,jj,jk) = ta(ji,jj,jk) + zta       ! add qsr trend to the temperature trend
139   
140#      if defined key_trdtra || defined key_trdmld
141                  ttrd(ji,jj,jk,7) = zta                  ! save the qsr trend
142#      endif
143               END DO
144            END DO
145            !                                             ! ===============
146         END DO                                           !   End of slab
147         !                                                ! ===============
148      ENDIF
149      !                                                ! =================== !
150      IF( lk_zco ) THEN                                !     z-coordinate    !
151         !                                             ! =================== !
152         !
153         !                                                ! ===============
154         DO jk = 1, nksr                                  ! Horizontal slab
155            !                                             ! ===============
156            zc0 = 1. / fse3t(1,1,jk)
157            DO jj = 2, jpjm1
158               DO ji = fs_2, fs_jpim1   ! vector opt.
159                  !                                       ! compute qsr forcing trend
160                  zta = qsr(ji,jj) * zc0 * ( gdsr(jk)*tmask(ji,jj,jk) - gdsr(jk+1)*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#      if defined key_trdtra || defined key_trdmld
165                  ttrd(ji,jj,jk,7) = zta                  ! save the qsr forcing trend
166#      endif
167               END DO
168            END DO
169            !                                             ! ===============
170         END DO                                           !   End of slab
171         !                                                ! ===============
172      ENDIF
173
174      IF( l_ctl .AND. lwp ) THEN         ! print mean trends (used for debugging)
175         zta = SUM( ta(2:jpim1,2:jpjm1,1:jpkm1) * tmask(2:jpim1,2:jpjm1,1:jpkm1) )
176         zta = SUM( ta * tmask ) 
177         WRITE(numout,*) ' qsr  - Ta: ', zta-t_ctl
178         t_ctl = zta 
179      ENDIF
180
181
182      ! 2. Substract qsr from the total heat flux q
183      ! -------------------------------------------
184      DO jj = 2, jpjm1
185         DO ji = fs_2, fs_jpim1   ! vector opt.
186            q(ji,jj) = qt(ji,jj) - qsr(ji,jj)
187         END DO
188      END DO
189
190   END SUBROUTINE tra_qsr
191
192
193   SUBROUTINE tra_qsr_init
194      !!----------------------------------------------------------------------
195      !!                  ***  ROUTINE tra_qsr_init  ***
196      !!
197      !! ** Purpose :   Initialization for the penetrative solar radiation
198      !!
199      !! ** Method  :   The profile of solar radiation within the ocean is set
200      !!      from two length scale of penetration (xsr1,xsr2) and a ratio
201      !!      (rabs). These parameters are read in the namqsr namelist. The
202      !!      default values correspond to clear water (type I in Jerlov'
203      !!      (1968) classification.
204      !!         called by tra_qsr at the first timestep (nit000)
205      !!
206      !! ** Action  : - initialize xsr1, xsr2 and rabs
207      !!
208      !! Reference :
209      !!   Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
210      !!
211      !! History :
212      !!   8.5  !  02-06  (G. Madec) Original code
213      !!----------------------------------------------------------------------
214      !! * Local declarations
215      INTEGER ::    jk,     &  ! dummy loop index
216                    indic      ! temporary integer
217      REAL(wp) ::   zdp1       ! temporary scalar
218
219      NAMELIST/namqsr/ ln_traqsr, rabs, xsi1, xsi2
220      !!----------------------------------------------------------------------
221
222      ! Read Namelist namqsr : ratio and length of penetration
223      ! --------------------
224      REWIND ( numnam )
225      READ   ( numnam, namqsr )
226
227      ! Parameter control and print
228      ! ---------------------------
229      IF( ln_traqsr  ) THEN
230        IF ( lwp ) THEN
231         WRITE(numout,*)
232         WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation'
233         WRITE(numout,*) '~~~~~~~~~~~~'
234         WRITE(numout,*) '          Namelist namqsr : set the parameter of penetration'
235         WRITE(numout,*) '             fraction associated with xsi   rabs   = ',rabs
236         WRITE(numout,*) '             first depth of extinction      xsi1   = ',xsi1
237         WRITE(numout,*) '             second depth of extinction     xsi2   = ',xsi2
238         WRITE(numout,*)
239        END IF
240      ELSE
241        IF ( lwp ) THEN
242         WRITE(numout,*)
243         WRITE(numout,*) 'tra_qsr_init : NO solar flux penetration'
244         WRITE(numout,*) '~~~~~~~~~~~~'
245        END IF
246      ENDIF
247
248      IF( rabs > 1.e0 .OR. rabs < 0.e0 .OR. xsi1 < 0.e0 .OR. xsi2 < 0.e0 ) THEN
249         IF(lwp) WRITE(numout,cform_err)
250         IF(lwp) WRITE(numout,*) '             0<rabs<1, 0<xsi1, or 0<xsi2 not satisfied'
251         nstop = nstop + 1
252      ENDIF
253
254
255      ! Initialization
256      ! --------------
257      IF( .NOT. lk_sco ) THEN
258         ! z-coordinate with or without partial step : same before last ocean w-level everywhere
259         DO jk = 1, jpk
260            zdp1 = -fsdepw(1,1,jk)
261            gdsr(jk) = ro0cpr * (  rabs  * EXP( zdp1/xsi1 ) + (1.-rabs) * EXP( zdp1/xsi2 )  )
262         END DO
263         indic = 0
264         DO jk = 1, jpk
265            IF( gdsr(jk) <= 1.e-15 .AND. indic == 0 ) THEN
266               gdsr(jk) = 0.e0
267               nksr = jk
268               !!bug Edmee chg res   nksr = jk - 1
269               indic = 1
270            ENDIF
271         END DO
272         nksr = MIN( nksr, jpkm1 )
273         IF(lwp) THEN
274            WRITE(numout,*)
275            WRITE(numout,*) '        - z-coordinate, level max of computation =', nksr
276            WRITE(numout,*) '             profile of coef. of penetration:'
277            WRITE(numout,"('              ',7e11.2)") ( gdsr(jk), jk = 1, nksr )
278            WRITE(numout,*)
279         ENDIF
280      ENDIF
281
282   END SUBROUTINE tra_qsr_init
283
284   !!======================================================================
285END MODULE traqsr
Note: See TracBrowser for help on using the repository browser.