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

Last change on this file since 145 was 106, checked in by opalod, 20 years ago

CT : UPDATE067 : Add control indices nictl, njctl used in SUM function output to compare mono versus multi procs runs
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.3 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) THEN         ! print mean trends (used for debugging)
175         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
176         WRITE(numout,*) ' qsr  - Ta: ', zta-t_ctl
177         t_ctl = zta 
178      ENDIF
179
180
181      ! 2. Substract qsr from the total heat flux q
182      ! -------------------------------------------
183      DO jj = 2, jpjm1
184         DO ji = fs_2, fs_jpim1   ! vector opt.
185            q(ji,jj) = qt(ji,jj) - qsr(ji,jj)
186         END DO
187      END DO
188
189   END SUBROUTINE tra_qsr
190
191
192   SUBROUTINE tra_qsr_init
193      !!----------------------------------------------------------------------
194      !!                  ***  ROUTINE tra_qsr_init  ***
195      !!
196      !! ** Purpose :   Initialization for the penetrative solar radiation
197      !!
198      !! ** Method  :   The profile of solar radiation within the ocean is set
199      !!      from two length scale of penetration (xsr1,xsr2) and a ratio
200      !!      (rabs). These parameters are read in the namqsr namelist. The
201      !!      default values correspond to clear water (type I in Jerlov'
202      !!      (1968) classification.
203      !!         called by tra_qsr at the first timestep (nit000)
204      !!
205      !! ** Action  : - initialize xsr1, xsr2 and rabs
206      !!
207      !! Reference :
208      !!   Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
209      !!
210      !! History :
211      !!   8.5  !  02-06  (G. Madec) Original code
212      !!----------------------------------------------------------------------
213      !! * Local declarations
214      INTEGER ::    jk,     &  ! dummy loop index
215                    indic      ! temporary integer
216      REAL(wp) ::   zdp1       ! temporary scalar
217
218      NAMELIST/namqsr/ ln_traqsr, rabs, xsi1, xsi2
219      !!----------------------------------------------------------------------
220
221      ! Read Namelist namqsr : ratio and length of penetration
222      ! --------------------
223      REWIND ( numnam )
224      READ   ( numnam, namqsr )
225
226      ! Parameter control and print
227      ! ---------------------------
228      IF( ln_traqsr  ) THEN
229        IF ( lwp ) THEN
230         WRITE(numout,*)
231         WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation'
232         WRITE(numout,*) '~~~~~~~~~~~~'
233         WRITE(numout,*) '          Namelist namqsr : set the parameter of penetration'
234         WRITE(numout,*) '             fraction associated with xsi   rabs   = ',rabs
235         WRITE(numout,*) '             first depth of extinction      xsi1   = ',xsi1
236         WRITE(numout,*) '             second depth of extinction     xsi2   = ',xsi2
237         WRITE(numout,*)
238        END IF
239      ELSE
240        IF ( lwp ) THEN
241         WRITE(numout,*)
242         WRITE(numout,*) 'tra_qsr_init : NO solar flux penetration'
243         WRITE(numout,*) '~~~~~~~~~~~~'
244        END IF
245      ENDIF
246
247      IF( rabs > 1.e0 .OR. rabs < 0.e0 .OR. xsi1 < 0.e0 .OR. xsi2 < 0.e0 ) THEN
248         IF(lwp) WRITE(numout,cform_err)
249         IF(lwp) WRITE(numout,*) '             0<rabs<1, 0<xsi1, or 0<xsi2 not satisfied'
250         nstop = nstop + 1
251      ENDIF
252
253
254      ! Initialization
255      ! --------------
256      IF( .NOT. lk_sco ) THEN
257         ! z-coordinate with or without partial step : same before last ocean w-level everywhere
258         DO jk = 1, jpk
259            zdp1 = -fsdepw(1,1,jk)
260            gdsr(jk) = ro0cpr * (  rabs  * EXP( zdp1/xsi1 ) + (1.-rabs) * EXP( zdp1/xsi2 )  )
261         END DO
262         indic = 0
263         DO jk = 1, jpk
264            IF( gdsr(jk) <= 1.e-15 .AND. indic == 0 ) THEN
265               gdsr(jk) = 0.e0
266               nksr = jk
267               !!bug Edmee chg res   nksr = jk - 1
268               indic = 1
269            ENDIF
270         END DO
271         nksr = MIN( nksr, jpkm1 )
272         IF(lwp) THEN
273            WRITE(numout,*)
274            WRITE(numout,*) '        - z-coordinate, level max of computation =', nksr
275            WRITE(numout,*) '             profile of coef. of penetration:'
276            WRITE(numout,"('              ',7e11.2)") ( gdsr(jk), jk = 1, nksr )
277            WRITE(numout,*)
278         ENDIF
279      ENDIF
280
281   END SUBROUTINE tra_qsr_init
282
283   !!======================================================================
284END MODULE traqsr
Note: See TracBrowser for help on using the repository browser.