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_tam.F90 in branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/TRA – NEMO

source: branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/TRA/traqsr_tam.F90 @ 3640

Last change on this file since 3640 was 3640, checked in by pabouttier, 10 years ago

Missing allocation/deallocation in TAM routines - See ticket #1013

  • Property svn:executable set to *
File size: 37.1 KB
Line 
1MODULE traqsr_tam
2#ifdef key_tam
3   !!======================================================================
4   !!                       ***  MODULE  traqsr_tam  ***
5   !! Ocean physics: solar radiation penetration in the top ocean levels
6   !!                Tangent and Adjoint Module
7   !!======================================================================
8   !! History :  OPA  !  1990-10  (B. Blanke)  Original code
9   !!            7.0  !  1991-11  (G. Madec)
10   !!                 !  1996-01  (G. Madec)  s-coordinates
11   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module
12   !!             -   !  2005-11  (G. Madec) zco, zps, sco coordinate
13   !!            3.2  !  2009-04  (G. Madec & NEMO team)
14   !! History of the TAM:
15   !!                 !  2008-05  (A. Vidard) Skeleton
16   !!            3.0  !  2008-09  (A. Vidard)   TAM of the 2005-11 version
17   !!            3.2  !  2010-03  (F. Vigilant) TAM of the 2009-11 version
18   !!            3.4  !  2012-07  (P.-A. Bouttier) 3.4 version
19   !!----------------------------------------------------------------------
20
21   !!----------------------------------------------------------------------
22   !!   tra_qsr      : trend due to the solar radiation penetration
23   !!   tra_qsr_init : solar radiation penetration initialization
24   !!----------------------------------------------------------------------
25   USE par_kind
26   USE par_oce
27   USE oce_tam
28   USE dom_oce
29   USE in_out_manager
30   USE fldread
31   USE sbc_oce
32   USE sbc_oce_tam
33   USE phycst
34   USE prtctl
35   USE gridrandom
36   USE dotprodfld
37   USE traqsr
38   USE trc_oce
39   USE trc_oce_tam
40   USE tstool_tam
41   USE lib_mpp
42   USE wrk_nemo
43   USE timing
44   USE restart
45   USE fldread
46   USE iom
47
48   IMPLICIT NONE
49   PRIVATE
50
51   PUBLIC   tra_qsr_tan      ! routine called by step_tam.F90 (ln_traqsr=T)
52   PUBLIC   tra_qsr_adj      ! routine called by step_tam.F90 (ln_traqsr=T)
53   PUBLIC   tra_qsr_init_tam
54   PUBLIC   tra_qsr_adj_tst  ! routine called by tst.F90
55
56   REAL(wp) :: xsi0r
57   REAL(wp) :: xsi1r
58   REAL(wp), DIMENSION(3,61) :: rkrgb
59   !! * Substitutions
60#  include "domzgr_substitute.h90"
61#  include "vectopt_loop_substitute.h90"
62
63CONTAINS
64
65   SUBROUTINE tra_qsr_tan( kt )
66      !!----------------------------------------------------------------------
67      !!                  ***  ROUTINE tra_qsr_tan  ***
68      !!
69      !! ** Purpose :   Compute the temperature trend due to the solar radiation
70      !!      penetration and add it to the general temperature trend.
71      !!
72      !! ** Method  : The profile of the solar radiation within the ocean is defined
73      !!      through 2 wavebands (rn_si0,rn_si1) or 3 wavebands (RGB) and a ratio rn_abs
74      !!      Considering the 2 wavebands case:
75      !!         I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) )
76      !!         The temperature trend associated with the solar radiation penetration
77      !!         is given by : zta = 1/e3t dk[ I ] / (rau0*Cp)
78      !!         At the bottom, boudary condition for the radiation is no flux :
79      !!      all heat which has not been absorbed in the above levels is put
80      !!      in the last ocean level.
81      !!         In z-coordinate case, the computation is only done down to the
82      !!      level where I(k) < 1.e-15 W/m2. In addition, the coefficients
83      !!      used for the computation are calculated one for once as they
84      !!      depends on k only.
85      !!
86      !! ** Action  : - update ta with the penetrative solar radiation trend
87      !!
88      !! Reference  : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
89      !!              Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516.
90      !!----------------------------------------------------------------------
91      INTEGER, INTENT(in) ::   kt        ! ocean time-step
92      !
93      !!
94      INTEGER  ::    ji, jj, jk          ! dummy loop indexes
95      INTEGER  ::   irgb                 ! temporary integers
96      REAL(wp) ::   zchl, zcoef, zfact, z1_e3t   ! temporary scalars
97      REAL(wp) ::   zc0, zc1, zc2, zc3   !    -         -
98      REAL(wp), POINTER, DIMENSION(:,:)     ::   zekb, zekg, zekr                      ! 2D workspace
99      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   ze0tl, ze1tl , ze2tl, ze3tl, zeatl    ! 3D workspace
100      !!----------------------------------------------------------------------
101      !
102      IF( nn_timing == 1 )  CALL timing_start('tra_qsr_tan')
103      !
104      CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        )
105      CALL wrk_alloc( jpi, jpj, jpk, ze0tl, ze1tl, ze2tl, ze3tl, zeatl )
106      !
107
108      IF( kt == nit000 ) THEN
109         IF(lwp) WRITE(numout,*)
110         IF(lwp) WRITE(numout,*) 'tra_qsr_tan : penetration of the surface solar radiation'
111         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
112         IF( .NOT.ln_traqsr )   RETURN
113      ENDIF
114      !                                        Set before qsr tracer content field
115      !                                        ***********************************
116      IF( kt == nit000 ) THEN                     ! Set the forcing field at nit000 - 1
117         !                                        ! -----------------------------------
118         IF( ln_rstart ) THEN !.AND.    &                    ! Restart: read in restart file
119              !& iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN
120            !IF(lwp) WRITE(numout,*) '          nit000-1 qsr tracer content forcing field red in the restart file'
121            zfact = 0.5e0
122            !CALL iom_get( numror, jpdom_autoglo, 'qsr_hc_b', qsr_hc_b_tl )   ! before heat content trend due to Qsr flux
123         ELSE                                           ! No restart or restart not found: Euler forward time stepping
124            zfact = 1.e0
125         ENDIF
126         qsr_hc_b_tl(:,:,:) = 0.e0
127      ELSE                                        ! Swap of forcing field
128         !                                        ! ---------------------
129         zfact = 0.5e0
130         qsr_hc_b_tl(:,:,:) = qsr_hc_tl(:,:,:)
131      ENDIF
132      !                                        Compute now qsr tracer content field
133      !
134      !                                           ! ============================================== !
135      IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN      !  bio-model fluxes  : all vertical coordinates  !
136         !                                        ! ============================================== !
137         DO jk = 1, jpkm1
138            qsr_hc_tl(:,:,jk) = ro0cpr * ( etot3_tl(:,:,jk) - etot3_tl(:,:,jk+1) )
139         END DO
140
141         DO jk = 1, jpkm1
142            DO jj = 2, jpjm1
143               DO ji = fs_2, fs_jpim1   ! vector opt.
144                  z1_e3t = zfact / fse3t(ji,jj,jk)
145                  tsa_tl(ji,jj,jk,jp_tem) = tsa_tl(ji,jj,jk,jp_tem) + ( qsr_hc_b_tl(ji,jj,jk) + qsr_hc_tl(ji,jj,jk) ) * z1_e3t
146               END DO
147            END DO
148         END DO
149         !                                        ! ============================================== !
150      ELSE                                        !  Ocean alone :
151         !                                        ! ============================================== !
152         !
153         !                                                ! ------------------------- !
154         IF( ln_qsr_rgb) THEN                             !  R-G-B  light penetration !
155            !                                             ! ------------------------- !
156            IF( nn_chldta == 1 .OR. lk_vvl ) THEN            !*  Variable Chlorophyll or ocean volume
157               !!! Set chlorophyl concentration
158               !!IF( nn_chldta ==1 ) THEN                             !*  Variable Chlorophyll
159                  !!!
160                  !!CALL fld_read( kt, 1, sf_chl )                         ! Read Chl data and provides it at the current time step
161                  !!!
162   !!!CDIR COLLAPSE
163   !!!CDIR NOVERRCHK
164                  !!DO jj = 1, jpj                                         ! Separation in R-G-B depending of the surface Chl
165   !!!CDIR NOVERRCHK
166                     !!DO ji = 1, jpi
167                        !!zchl = MIN( 10.0_wp , MAX( 0.03_wp, sf_chl(1)%fnow(ji,jj) ) )
168                        !!irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
169                        !!zekb(ji,jj) = rkrgb(1,irgb)
170                        !!zekg(ji,jj) = rkrgb(2,irgb)
171                        !!zekr(ji,jj) = rkrgb(3,irgb)
172                     !!END DO
173                  !!END DO
174               !ELSE
175                  !zchl = 0.05                                     ! constant chlorophyll
176                  !irgb = NINT( 41 + 20.*LOG10( zchl ) + 1.e-15 )
177                  !zekb(:,:) = rkrgb(1,irgb)                       ! Separation in R-G-B depending of the chlorophyll
178                  !zekg(:,:) = rkrgb(2,irgb)
179                  !zekr(:,:) = rkrgb(3,irgb)
180               !ENDIF
181               !
182               !zcoef  = ( 1.0_wp - rn_abs ) / 3.0_wp                        ! equi-partition in R-G-B
183               !ze0tl(:,:,1) = rn_abs  * qsr_tl(:,:)
184               !ze1tl(:,:,1) = zcoef * qsr_tl(:,:)
185               !ze2tl(:,:,1) = zcoef * qsr_tl(:,:)
186               !ze3tl(:,:,1) = zcoef * qsr_tl(:,:)
187               !zeatl(:,:,1) =         qsr_tl(:,:)
188               !!
189               !DO jk = 2, nksr+1
190!!CDIR NOVERRCHK
191                  !DO jj = 1, jpj
192!!CDIR NOVERRCHK
193                     !DO ji = 1, jpi
194                        !zc0tl = ze0tl(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r     )
195                        !zc1tl = ze1tl(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekb(ji,jj) )
196                        !zc2tl = ze2tl(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekg(ji,jj) )
197                        !zc3tl = ze3tl(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekr(ji,jj) )
198                        !ze0tl(ji,jj,jk) = zc0tl
199                        !ze1tl(ji,jj,jk) = zc1tl
200                        !ze2tl(ji,jj,jk) = zc2tl
201                        !ze3tl(ji,jj,jk) = zc3tl
202                        !zeatl(ji,jj,jk) = ( zc0tl + zc1tl + zc2tl + zc3tl ) * tmask(ji,jj,jk)
203                     !END DO
204                  !END DO
205               !END DO
206               !!
207               !DO jk = 1, nksr                                        ! compute and add qsr trend to ta
208                  !qsr_tl(:,:) = ro0cpr * ( zeatl(:,:,jk) - zeatl(:,:,jk+1) )
209               !END DO
210               !zeatl(:,:,nksr+1:jpk) = 0.0_wp     ! below 400m set to zero
211               !!
212               CALL ctl_stop('tra_qsr_tan: key_vvl or non-constant chlorophyll management(nn_chldta = 1) &
213                          &   not implemented in TAM yet')
214            ELSE                                                 !*  Constant Chlorophyll
215               DO jk = 1, nksr
216                  qsr_hc_tl(:,:,jk) = etot3_tl(:,:,jk) * qsr(:,:) + etot3(:,:,jk) * qsr_tl(:,:)
217               END DO
218            ENDIF
219
220         ENDIF
221         !                                                ! ------------------------- !
222         IF( ln_qsr_2bd ) THEN                            !  2 band light penetration !
223            !                                             ! ------------------------- !
224            !
225            IF( lk_vvl ) THEN                                  !* variable volume
226               !zz0   =        rn_abs   * ro0cpr
227               !zz1   = ( 1. - rn_abs ) * ro0cpr
228               !DO jk = 1, nksr                    ! solar heat absorbed at T-point in the top 400m
229                  !DO jj = 1, jpj
230                     !DO ji = 1, jpi
231                        !zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r )
232                        !zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
233                        !qsr_hc_tl(ji,jj,jk) = qsr_tl(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) )
234                     !END DO
235                  !END DO
236               !END DO
237               CALL ctl_stop('tra_qsr_tan: key_vvl or chlorophyll management not implemented in TAM yet')
238            ELSE                                               !* constant volume: coef. computed one for all
239               DO jk = 1, nksr
240                  DO jj = 2, jpjm1
241                     DO ji = fs_2, fs_jpim1   ! vector opt.
242                        qsr_hc_tl(ji,jj,jk) =  etot3_tl(ji,jj,jk) * qsr(ji,jj) + etot3(ji,jj,jk) * qsr_tl(ji,jj)
243                     END DO
244                  END DO
245               END DO
246               !
247            ENDIF
248            !
249         ENDIF
250         DO jk = 1, nksr
251            DO jj = 2, jpjm1
252               DO ji = fs_2, fs_jpim1   ! vector opt.
253                  z1_e3t = zfact / fse3t(ji,jj,jk)
254                  tsa_tl(ji,jj,jk,jp_tem) = tsa_tl(ji,jj,jk,jp_tem) + ( qsr_hc_b_tl(ji,jj,jk) + qsr_hc_tl(ji,jj,jk) ) * z1_e3t
255               END DO
256            END DO
257         END DO
258         !
259      ENDIF
260      !
261      CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        )
262      CALL wrk_dealloc( jpi, jpj, jpk, ze0tl, ze1tl, ze2tl, ze3tl, zeatl )
263      !
264      IF( nn_timing == 1 )  CALL timing_stop('tra_qsr_tan')
265      !
266   END SUBROUTINE tra_qsr_tan
267   SUBROUTINE tra_qsr_adj( kt )
268      !!----------------------------------------------------------------------
269      !!                  ***  ROUTINE tra_qsr_adj  ***
270      !!
271      !! ** Purpose :   Compute the temperature trend due to the solar radiation
272      !!      penetration and add it to the general temperature trend.
273      !!
274      !! ** Method  : The profile of the solar radiation within the ocean is defined
275      !!      through 2 wavebands (rn_si0,rn_si1) or 3 wavebands (RGB) and a ratio rn_abs
276      !!      Considering the 2 wavebands case:
277      !!         I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) )
278      !!         The temperature trend associated with the solar radiation penetration
279      !!         is given by : zta = 1/e3t dk[ I ] / (rau0*Cp)
280      !!         At the bottom, boudary condition for the radiation is no flux :
281      !!      all heat which has not been absorbed in the above levels is put
282      !!      in the last ocean level.
283      !!         In z-coordinate case, the computation is only done down to the
284      !!      level where I(k) < 1.e-15 W/m2. In addition, the coefficients
285      !!      used for the computation are calculated one for once as they
286      !!      depends on k only.
287      !!
288      !! ** Action  : - update ta with the penetrative solar radiation trend
289      !!
290      !! Reference  : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
291      !!              Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516.
292      !!----------------------------------------------------------------------
293      !!
294      INTEGER, INTENT(in) ::   kt     ! ocean time-step
295      !
296      !!
297      INTEGER  ::    ji, jj, jk       ! dummy loop indexes
298      INTEGER  ::   irgb                 ! temporary integers
299      REAL(wp) ::   zchl, zcoef, zfact, z1_e3t   ! temporary scalars
300      REAL(wp) ::   zc0, zc1, zc2, zc3, zz0, zz1   !    -         -
301      REAL(wp), POINTER, DIMENSION(:,:)     ::   zekb, zekg, zekr                      ! 2D workspace
302      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ze0ad, ze1ad , ze2ad, ze3ad, zeaad    ! 3D workspace
303      !!----------------------------------------------------------------------
304      !
305      IF( nn_timing == 1 )  CALL timing_start('tra_qsr_adj')
306      !
307      CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        )
308      CALL wrk_alloc( jpi, jpj, jpk, ze0ad, ze1ad, ze2ad, ze3ad, zeaad )
309      !
310      IF( kt == nitend ) THEN
311         IF(lwp) WRITE(numout,*)
312         IF(lwp) WRITE(numout,*) 'tra_qsr_adj : penetration of the surface solar radiation'
313         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
314         IF( .NOT.ln_traqsr )   RETURN
315      ENDIF
316      !                                        Set before qsr tracer content field
317      !                                        ***********************************
318      IF( kt == nit000 ) THEN                     ! Set the forcing field at nit000 - 1
319         !                                        ! -----------------------------------
320         IF( ln_rstart ) THEN !.AND.    &                    ! Restart: read in restart file
321              !& iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN
322            !IF(lwp) WRITE(numout,*) '          nit000-1 qsr tracer content forcing field red in the restart file'
323            zfact = 0.5e0
324            !CALL iom_get( numror, jpdom_autoglo, 'qsr_hc_b', qsr_hc_b_ad )   ! before heat content trend due to Qsr flux
325         ELSE                                           ! No restart or restart not found: Euler forward time stepping
326            zfact = 1.e0
327         ENDIF
328      ELSE                                        ! Swap of forcing field
329         !                                        ! ---------------------
330         zfact = 0.5e0
331      ENDIF
332      !                                           ! ============================================== !
333      IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN      !  bio-model fluxes  : all vertical coordinates  !
334         !                                        ! ============================================== !
335         DO jk = jpkm1, 1, -1
336            DO jj = 2, jpjm1
337               DO ji = fs_2, fs_jpim1   ! vector opt.
338                  z1_e3t = zfact / fse3t(ji,jj,jk)
339                  qsr_hc_b_ad(ji,jj,jk) = qsr_hc_b_ad(ji,jj,jk) + tsa_ad(ji,jj,jk,jp_tem) * z1_e3t
340                  qsr_hc_ad(ji,jj,jk) = qsr_hc_ad(ji,jj,jk) + tsa_ad(ji,jj,jk,jp_tem) * z1_e3t
341               END DO
342            END DO
343         END DO
344         DO jk = jpkm1, 1, -1
345            etot3_ad(:,:,jk)   = etot3_ad(:,:,jk) + ro0cpr * qsr_hc_ad(:,:,jk)
346            etot3_ad(:,:,jk+1) = etot3_ad(:,:,jk+1) - ro0cpr * qsr_hc_ad(:,:,jk)
347         END DO
348         !                                        ! ============================================== !
349      ELSE                                        !  Ocean alone :
350         !                                        ! ============================================== !
351         !
352         DO jk = 1, nksr
353            DO jj = 2, jpjm1
354               DO ji = fs_2, fs_jpim1   ! vector opt.
355                  z1_e3t = zfact / fse3t(ji,jj,jk)
356                  qsr_hc_b_ad(ji,jj,jk) = qsr_hc_b_ad(ji,jj,jk) + tsa_ad(ji,jj,jk,jp_tem) * z1_e3t
357                  qsr_hc_ad(ji,jj,jk)   = qsr_hc_ad(ji,jj,jk)   + tsa_ad(ji,jj,jk,jp_tem) * z1_e3t
358               END DO
359            END DO
360         END DO
361         !                                                ! ------------------------- !
362         IF( ln_qsr_2bd ) THEN                            !  2 band light penetration !
363            !                                             ! ------------------------- !
364            !
365            IF( lk_vvl ) THEN                                  !* variable volume
366               !zz0   =        rn_abs   * ro0cpr
367               !zz1   = ( 1. - rn_abs ) * ro0cpr
368               !DO jk = nksr, 1, -1                    ! solar heat absorbed at T-point in the top 400m
369                  !DO jj = 1, jpj
370                     !DO ji = 1, jpi
371                        !zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r )
372                        !zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
373                        !qsr_ad(ji,jj) = qsr_hc_ad(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) )
374                     !END DO
375                  !END DO
376               !END DO
377               CALL ctl_stop('tra_qsr_adj: key_vvl or chlorophyll management not implemented in TAM yet')
378            ELSE                                               !* constant volume: coef. computed one for all
379               DO jk = 1, nksr
380                  DO jj = 2, jpjm1
381                     DO ji = fs_2, fs_jpim1   ! vector opt.
382                        etot3_ad(ji,jj,jk) = etot3_ad(ji,jj,jk) + qsr(ji,jj) * qsr_hc_ad(ji,jj,jk)
383                        qsr_ad(ji,jj)      = qsr_ad(ji,jj) + etot3(ji,jj,jk) * qsr_hc_ad(ji,jj,jk)
384                        qsr_hc_ad(ji,jj,jk) = 0._wp
385                     END DO
386                  END DO
387               END DO
388               !
389            ENDIF
390            !
391         ENDIF
392         !
393         !                                                ! ------------------------- !
394         IF( ln_qsr_rgb) THEN                             !  R-G-B  light penetration !
395            !                                             ! ------------------------- !
396            ! Set chlorophyl concentration
397            IF( nn_chldta == 1 .OR. lk_vvl ) THEN            !*  Variable Chlorophyll or ocean volume
398               !!!
399               !!IF( nn_chldta ==1 ) THEN                             !*  Variable Chlorophyll
400                  !!zc0ad = 0.0_wp; zc1ad = 0.0_wp; zc2ad = 0.0_wp; zc3ad = 0.0_wp
401                  !!ze0ad(:,:,:) = 0.0_wp; ze1ad(:,:,:) = 0.0_wp; ze2ad(:,:,:) = 0.0_wp; ze3ad(:,:,:) = 0.0_wp
402                  !!zeaad(:,:,:) = 0.0_wp
403                  !!!
404                  !!CALL fld_read( kt, 1, sf_chl )                         ! Read Chl data and provides it at the current time step
405                  !!!
406   !!!CDIR COLLAPSE
407   !!!CDIR NOVERRCHK
408                  !!DO jj = 1, jpj                                         ! Separation in R-G-B depending of the surface Chl
409   !!!CDIR NOVERRCHK
410                     !!DO ji = 1, jpi
411                        !!zchl = MIN( 10.0_wp , MAX( 0.03_wp, sf_chl(1)%fnow(ji,jj) ) )
412                        !!irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
413                        !!zekb(ji,jj) = rkrgb(1,irgb)
414                        !!zekg(ji,jj) = rkrgb(2,irgb)
415                        !!zekr(ji,jj) = rkrgb(3,irgb)
416                     !!END DO
417                  !!END DO
418               !ELSE
419                  !zchl = 0.05                                     ! constant chlorophyll
420                  !irgb = NINT( 41 + 20.*LOG10( zchl ) + 1.e-15 )
421                  !zekb(:,:) = rkrgb(1,irgb)                       ! Separation in R-G-B depending of the chlorophyll
422                  !zekg(:,:) = rkrgb(2,irgb)
423                  !zekr(:,:) = rkrgb(3,irgb)
424               !ENDIF
425               !!
426               !zcoef  = ( 1.0_wp - rn_abs ) / 3.0_wp
427
428               !zeaad(:,:,nksr+1:jpk) = 0.0_wp     ! below 400m set to zero
429               !!
430               !DO jk = 1, nksr                                        ! compute and add qsr trend to ta
431                  !zeaad(:,:,jk  ) =   ro0cpr * qsr_hc_ad(:,:,jk)
432                  !zeaad(:,:,jk+1) = - ro0cpr * qsr_hc_ad(:,:,jk)
433               !END DO
434               !!
435               !DO jk = nksr+1, 2, -1
436!!CDIR NOVERRCHK
437                  !DO jj = 1, jpj
438!!CDIR NOVERRCHK
439                     !DO ji = 1, jpi
440                        !zc0ad = zc0ad + zeaad(ji,jj,jk) * tmask(ji,jj,jk)
441                        !zc1ad = zc1ad + zeaad(ji,jj,jk) * tmask(ji,jj,jk)
442                        !zc2ad = zc2ad + zeaad(ji,jj,jk) * tmask(ji,jj,jk)
443                        !zc3ad = zc3ad + zeaad(ji,jj,jk) * tmask(ji,jj,jk)
444                        !zeaad(ji,jj,jk) = 0.0_wp
445                        !zc0ad = zc0ad + ze0ad(ji,jj,jk)
446                        !zc1ad = zc1ad + ze1ad(ji,jj,jk)
447                        !zc2ad = zc2ad + ze2ad(ji,jj,jk)
448                        !zc3ad = zc3ad + ze3ad(ji,jj,jk)
449                        !ze0ad(ji,jj,jk) = 0.0_wp
450                        !ze1ad(ji,jj,jk) = 0.0_wp
451                        !ze2ad(ji,jj,jk) = 0.0_wp
452                        !ze3ad(ji,jj,jk) = 0.0_wp
453                        !ze0ad(ji,jj,jk-1) = ze0ad(ji,jj,jk-1) + zc0ad * EXP( - fse3t(ji,jj,jk-1) * xsi0r     )
454                        !ze1ad(ji,jj,jk-1) = ze1ad(ji,jj,jk-1) + zc1ad * EXP( - fse3t(ji,jj,jk-1) * zekb(ji,jj) )
455                        !ze2ad(ji,jj,jk-1) = ze2ad(ji,jj,jk-1) + zc2ad * EXP( - fse3t(ji,jj,jk-1) * zekg(ji,jj) )
456                        !ze3ad(ji,jj,jk-1) = ze3ad(ji,jj,jk-1) + zc3ad * EXP( - fse3t(ji,jj,jk-1) * zekr(ji,jj) )
457                        !zc0ad = 0.0_wp
458                        !zc1ad = 0.0_wp
459                        !zc2ad = 0.0_wp
460                        !zc3ad = 0.0_wp
461                     !END DO
462                  !END DO
463               !END DO
464               !!
465               !qsr_ad(:,:) = qsr_ad(:,:) + zeaad(:,:,1)
466               !qsr_ad(:,:) = qsr_ad(:,:) + zcoef  * ze3ad(:,:,1)
467               !qsr_ad(:,:) = qsr_ad(:,:) + zcoef  * ze2ad(:,:,1)
468               !qsr_ad(:,:) = qsr_ad(:,:) + zcoef  * ze1ad(:,:,1)
469               !qsr_ad(:,:) = qsr_ad(:,:) + rn_abs * ze0ad(:,:,1)
470               !!
471               CALL ctl_stop('tra_qsr_adj: key_vvl or chlorophyll management not implemented in TAM yet')
472            ELSE                                                 !*  Constant Chlorophyll
473               DO jk = 1, nksr
474                  etot3_ad(:,:,jk) = etot3_ad(:,:,jk) + qsr_hc_ad(:,:,jk) * qsr(:,:)
475                  qsr_ad(  :,:   ) = qsr_ad(:,:)      + qsr_hc_ad(:,:,jk) * etot3(:,:,jk)
476                  qsr_hc_ad(:,:,jk) = 0._wp
477               END DO
478            ENDIF
479         ENDIF
480      ENDIF
481      IF ( kt /= nit000 ) THEN
482         qsr_hc_ad(:,:,:) = qsr_hc_ad(:,:,:) + qsr_hc_b_ad(:,:,:)
483      ENDIF
484      qsr_hc_b_ad(:,:,:) = 0._wp
485
486      CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        )
487      CALL wrk_dealloc( jpi, jpj, jpk, ze0ad, ze1ad, ze2ad, ze3ad, zeaad )
488
489      IF( nn_timing == 1 )  CALL timing_stop('tra_qsr_adj')
490
491   END SUBROUTINE tra_qsr_adj
492   SUBROUTINE tra_qsr_adj_tst ( kumadt )
493      !!-----------------------------------------------------------------------
494      !!
495      !!          ***  ROUTINE tra_sbc_adj_tst : TEST OF tra_sbc_adj  ***
496      !!
497      !! ** Purpose : Test the adjoint routine.
498      !!
499      !! ** Method  : Verify the scalar product
500      !!
501      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
502      !!
503      !!              where  L   = tangent routine
504      !!                     L^T = adjoint routine
505      !!                     W   = diagonal matrix of scale factors
506      !!                     dx  = input perturbation (random field)
507      !!                     dy  = L dx
508      !!
509      !! History :
510      !!        ! 08-08 (A. Vidard)
511      !!-----------------------------------------------------------------------
512      !! * Modules used
513
514      !! * Arguments
515      INTEGER, INTENT(IN) :: &
516         & kumadt             ! Output unit
517
518      INTEGER ::  &
519         & jstp,  &
520         & ji,    &        ! dummy loop indices
521         & jj,    &
522         & jk
523      !! * Local declarations
524      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
525         & zta_tlin,     &! Tangent input : after temperature
526         & zta_tlout,    &! Tangent output: after temperature
527         & zta_adout,    &! Adjoint output: after temperature
528         & zta_adin,     &! Adjoint input : after temperature
529         & zqsr_hc_tlin,     &! qsr_hcngent input : after temperature
530         & zqsr_hc_tlout,    &! qsr_hcngent output: after temperature
531         & zqsr_hc_adout,    &! Adjoint output: after temperature
532         & zqsr_hc_adin,     &! Adjoint input : after temperature
533         & zqsr_hc_b_tlout,    &! qsr_hc_bngent output: after temperature
534         & zqsr_hc_b_adin,     &! Adjoint input : after temperature
535         & zetot3_tlin,  &! Tangent input
536         & zetot3_adout, &! Adjoint output
537         & zta,          & ! temporary after temperature
538         & zqsr_hc,          & ! temporary after temperature
539         & zqsr_hc_b,          & ! temporary after temperature
540         & zetot3          ! temporary
541      REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
542         & zqsr_tlin,    &! Tangent input : solar radiation (w/m2)
543         & zqsr_adout,   &! Adjoint output: solar radiation (w/m2)
544         & zqsr           ! temporary solar radiation (w/m2)
545      REAL(KIND=wp) ::       &
546         & zsp1,             & ! scalar product involving the tangent routine
547         & zsp2,             & ! scalar product involving the adjoint routine
548         & zsp2_1,           & ! scalar product involving the adjoint routine
549         & zsp2_2,           & ! scalar product involving the adjoint routine
550         & zsp2_3              ! scalar product involving the adjoint routine
551      CHARACTER(LEN=14) :: &
552         & cl_name
553
554      ALLOCATE( &
555         & zta_tlin(jpi,jpj,jpk),     &
556         & zta_tlout(jpi,jpj,jpk),    &
557         & zta_adout(jpi,jpj,jpk),    &
558         & zta_adin(jpi,jpj,jpk),     &
559         & zta(jpi,jpj,jpk),          &
560         & zqsr_hc_tlin(jpi,jpj,jpk),     &
561         & zqsr_hc_tlout(jpi,jpj,jpk),    &
562         & zqsr_hc_adout(jpi,jpj,jpk),    &
563         & zqsr_hc_adin(jpi,jpj,jpk),     &
564         & zqsr_hc(jpi,jpj,jpk),          &
565         & zqsr_hc_b_tlout(jpi,jpj,jpk),    &
566         & zqsr_hc_b_adin(jpi,jpj,jpk),     &
567         & zqsr_hc_b(jpi,jpj,jpk),          &
568         & zqsr_tlin(jpi,jpj),        &
569         & zqsr_adout(jpi,jpj),       &
570         & zetot3_tlin(jpi,jpj,jpk),  &
571         & zetot3_adout(jpi,jpj,jpk), &
572         & zqsr(jpi,jpj),             &
573         & zetot3(jpi,jpj,jpk)        &
574         & )
575      ! Initialize the reference state
576      qsr(:,:) = 1.0_wp ! ???
577      !Initialize etot3 to non-zero value until traj(nit000-1) is fixed
578      etot3(:,:,1) = 2.e-8   ; etot3(:,:,2) = 1.5e-9; etot3(:,:,3) = 8.5e-10
579      etot3(:,:,4) = 5.4e-10 ; etot3(:,:,5) = 3.5e-10; etot3(:,:,6:jpk) = 0.0_wp
580      ! Initialize random field standard deviations
581      !=============================================================
582      ! 1) dx = ( T ) and dy = ( T )
583      !=============================================================
584
585      !--------------------------------------------------------------------
586      ! Reset the tangent and adjoint variables
587      !--------------------------------------------------------------------
588      zta_tlin(:,:,:)        = 0.0_wp
589      zta_tlout(:,:,:)       = 0.0_wp
590      zta_adout(:,:,:)       = 0.0_wp
591      zta_adin(:,:,:)        = 0.0_wp
592      zqsr_hc_tlin(:,:,:)    = 0.0_wp
593      zqsr_hc_tlout(:,:,:)   = 0.0_wp
594      zqsr_hc_adout(:,:,:)   = 0.0_wp
595      zqsr_hc_adin(:,:,:)    = 0.0_wp
596      zqsr_hc_b_tlout(:,:,:) = 0.0_wp
597      zqsr_hc_b_adin(:,:,:)  = 0.0_wp
598      zqsr_adout(:,:)        = 0.0_wp
599      zqsr_tlin(:,:)         = 0.0_wp
600      zetot3_tlin(:,:,:)     = 0.0_wp
601      zetot3_adout(:,:,:)    = 0.0_wp
602      tsa_ad(:,:,:,jp_tem)   = 0.0_wp
603      qsr_ad(:,:)            = 0.0_wp
604      qsr_hc_ad(:,:,:)       = 0.0_wp
605      qsr_hc_b_ad(:,:,:)     = 0.0_wp
606      etot3_ad(:,:,:)        = 0.0_wp
607
608      CALL grid_random(  zqsr   , 'T', 0.0_wp, stdqsr )
609      CALL grid_random(  zqsr_hc, 'T', 0.0_wp, stdqsr )
610      CALL grid_random(  zta    , 'T', 0.0_wp, stdt )
611      CALL grid_random(  zetot3 , 'T', 0.0_wp, stdt )
612      DO jk = 1, jpk
613         DO jj = nldj, nlej
614            DO ji = nldi, nlei
615               zta_tlin(ji,jj,jk) = zta(ji,jj,jk)
616            END DO
617         END DO
618      END DO
619      DO jk = 1, jpk
620         DO jj = nldj, nlej
621            DO ji = nldi, nlei
622               zqsr_hc_tlin(ji,jj,jk) = zqsr_hc(ji,jj,jk)
623            END DO
624         END DO
625      END DO
626      DO jk = 1, jpk
627         DO jj = nldj, nlej
628            DO ji = nldi, nlei
629               zetot3_tlin(ji,jj,jk) = zetot3(ji,jj,jk)
630            END DO
631         END DO
632      END DO
633      DO jj = nldj, nlej
634         DO ji = nldi, nlei
635            zqsr_tlin(ji,jj)  = zqsr(ji,jj)
636         END DO
637      END DO
638      ! Test for time steps nit000 and nit000 + 1 (the matrix changes)
639      DO jstp = nit000, nit000 + 1
640         !--------------------------------------------------------------------
641         ! Call the tangent routine: dy = L dx
642         !--------------------------------------------------------------------
643
644         tsa_tl(:,:,:,jp_tem) = zta_tlin(:,:,:)
645         etot3_tl(:,:,:)      = zetot3_tlin(:,:,:)
646         qsr_tl(:,:)          = zqsr_tlin(:,:)
647         qsr_hc_tl(:,:,:)     = zqsr_hc_tlin(:,:,:)
648
649         CALL tra_qsr_tan( jstp )
650
651         zta_tlout(:,:,:)       = tsa_tl(:,:,:,jp_tem)
652         zqsr_hc_tlout(:,:,:)   = qsr_hc_tl(:,:,:)
653         zqsr_hc_b_tlout(:,:,:) = qsr_hc_b_tl(:,:,:)
654
655         !--------------------------------------------------------------------
656         ! Initialize the adjoint variables: dy^* = W dy
657         !--------------------------------------------------------------------
658
659         DO jk = 1, jpk
660            DO jj = nldj, nlej
661               DO ji = nldi, nlei
662                  zta_adin(ji,jj,jk) = zta_tlout(ji,jj,jk) &
663                     &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
664                     &               * tmask(ji,jj,jk)
665               END DO
666            END DO
667         END DO
668         DO jk = 1, jpk
669            DO jj = nldj, nlej
670               DO ji = nldi, nlei
671                  zqsr_hc_adin(ji,jj,jk) = zqsr_hc_tlout(ji,jj,jk) &
672                     &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
673                     &               * tmask(ji,jj,jk)
674               END DO
675            END DO
676         END DO
677         DO jk = 1, jpk
678            DO jj = nldj, nlej
679               DO ji = nldi, nlei
680                  zqsr_hc_b_adin(ji,jj,jk) = zqsr_hc_b_tlout(ji,jj,jk) &
681                     &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
682                     &               * tmask(ji,jj,jk)
683               END DO
684            END DO
685         END DO
686
687         !--------------------------------------------------------------------
688         ! Compute the scalar product: ( L dx )^T W dy
689         !--------------------------------------------------------------------
690
691         zsp1 = DOT_PRODUCT( zta_tlout, zta_adin ) &
692            &  + DOT_PRODUCT( zqsr_hc_tlout, zqsr_hc_adin ) &
693            &  + DOT_PRODUCT( zqsr_hc_b_tlout, zqsr_hc_b_adin )
694
695         !--------------------------------------------------------------------
696         ! Call the adjoint routine: dx^* = L^T dy^*
697         !--------------------------------------------------------------------
698
699         etot3_ad(:,:,:)      = 0.0_wp
700         qsr_ad(:,:)          = 0.0_wp
701         tsa_ad(:,:,:,jp_tem) = zta_adin(:,:,:)
702         qsr_hc_ad(:,:,:)     = zqsr_hc_adin(:,:,:)
703         qsr_hc_b_ad(:,:,:)   = zqsr_hc_b_adin(:,:,:)
704
705         CALL tra_qsr_adj( jstp )
706
707         zta_adout(:,:,:)       = tsa_ad(:,:,:,jp_tem)
708         zetot3_adout(:,:,:)    = etot3_ad(:,:,:)
709         zqsr_adout(:,:)        = qsr_ad(:,:)
710         zqsr_hc_adout(:,:,:)   = qsr_hc_ad(:,:,:)
711
712         !--------------------------------------------------------------------
713         ! Compute the scalar product: dx^T L^T W dy
714         !--------------------------------------------------------------------
715
716         zsp2_1 = DOT_PRODUCT( zta_tlin    , zta_adout     )
717         zsp2_1 = zsp2_1 + DOT_PRODUCT( zqsr_hc_tlin    , zqsr_hc_adout   )
718         zsp2_2 = DOT_PRODUCT( zqsr_tlin   , zqsr_adout    )
719         zsp2_3 = DOT_PRODUCT( zetot3_tlin , zetot3_adout  )
720
721         zsp2 = zsp2_1 + zsp2_2 + zsp2_3
722
723         ! Compare the scalar products
724
725         ! 14 char:   '12345678901234'
726         IF (jstp == nit000) THEN
727            cl_name = 'tra_qsr_adj  1'
728         ELSE
729            cl_name = 'tra_qsr_adj  2'
730         END IF
731         CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
732      END DO
733
734      DEALLOCATE( &
735         & zta_tlin,        &
736         & zta_tlout,       &
737         & zta_adout,       &
738         & zta_adin,        &
739         & zta,             &
740         & zqsr_hc_tlin,    &
741         & zqsr_hc_tlout,   &
742         & zqsr_hc_adout,   &
743         & zqsr_hc_adin,    &
744         & zqsr_hc,         &
745         & zqsr_hc_b_tlout, &
746         & zqsr_hc_b_adin,  &
747         & zqsr_hc_b,       &
748         & zqsr_adout,      &
749         & zqsr_tlin,       &
750         & zqsr             &
751         & )
752
753      !
754   END SUBROUTINE tra_qsr_adj_tst
755   SUBROUTINE tra_qsr_init_tam
756      !!----------------------------------------------------------------------
757      !!                  ***  ROUTINE tra_qsr_init_tan  ***
758      !!
759      !! ** Purpose :   Initialization for the penetrative solar radiation
760      !!
761      !! ** Method  :   The profile of solar radiation within the ocean is set
762      !!      from two length scale of penetration (rn_si0,rn_si1) and a ratio
763      !!      (rn_abs). These parameters are read in the namtra_qsr namelist. The
764      !!      default values correspond to clear water (type I in Jerlov'
765      !!      (1968) classification.
766      !!         called by tra_qsr at the first timestep (nit000)
767      !!
768      !! ** Action  : - initialize rn_si0, rn_si1 and rn_abs
769      !!
770      !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
771      !!----------------------------------------------------------------------
772
773      IF( ln_traqsr  ) THEN      !  Initialisation of Light Penetration  !
774         !                       ! ===================================== !
775         !
776         !                                ! ---------------------------------- !
777         IF( ln_qsr_rgb ) THEN            !  Red-Green-Blue light penetration  !
778            !                             ! ---------------------------------- !
779            etot3_tl(:,:,:) = 0.0_wp
780            etot3_ad(:,:,:) = 0.0_wp
781            !
782         ENDIF
783            !                             ! ---------------------------------- !
784         IF( ln_qsr_2bd ) THEN            !    2 bands    light penetration    !
785            !                             ! ---------------------------------- !
786            etot3_tl(:,:,:) = 0.0_wp
787            etot3_ad(:,:,:) = 0.0_wp
788            !
789         ENDIF
790         !
791      ENDIF
792      !
793   END SUBROUTINE tra_qsr_init_tam
794
795   !!======================================================================
796#endif
797END MODULE traqsr_tam
Note: See TracBrowser for help on using the repository browser.