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

Context Navigation

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

Last change on this file since 3627 was 3627, checked in by rblod, 11 years ago
Correct long lines in TAM 4.3 see ticket #1010
Property svn:executable set to ``*
File size: 37.0 KB

Line
1	MODULE 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
63	CONTAINS
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_absEXP(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) * ( zc0tmask(ji,jj,jk) - zc1tmask(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_absEXP(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) * ( zc0tmask(ji,jj,jk) - zc1tmask(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	IF( nn_timing == 1 ) CALL timing_stop('tra_qsr_adj')
487
488	END SUBROUTINE tra_qsr_adj
489	SUBROUTINE tra_qsr_adj_tst ( kumadt )
490	!!-----------------------------------------------------------------------
491	!!
492	!! * ROUTINE tra_sbc_adj_tst : TEST OF tra_sbc_adj *
493	!!
494	!! ** Purpose : Test the adjoint routine.
495	!!
496	!! ** Method : Verify the scalar product
497	!!
498	!! ( L dx )^T W dy = dx^T L^T W dy
499	!!
500	!! where L = tangent routine
501	!! L^T = adjoint routine
502	!! W = diagonal matrix of scale factors
503	!! dx = input perturbation (random field)
504	!! dy = L dx
505	!!
506	!! History :
507	!! ! 08-08 (A. Vidard)
508	!!-----------------------------------------------------------------------
509	!! * Modules used
510
511	!! * Arguments
512	INTEGER, INTENT(IN) :: &
513	& kumadt ! Output unit
514
515	INTEGER :: &
516	& jstp, &
517	& ji, & ! dummy loop indices
518	& jj, &
519	& jk
520	!! * Local declarations
521	REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
522	& zta_tlin, &! Tangent input : after temperature
523	& zta_tlout, &! Tangent output: after temperature
524	& zta_adout, &! Adjoint output: after temperature
525	& zta_adin, &! Adjoint input : after temperature
526	& zqsr_hc_tlin, &! qsr_hcngent input : after temperature
527	& zqsr_hc_tlout, &! qsr_hcngent output: after temperature
528	& zqsr_hc_adout, &! Adjoint output: after temperature
529	& zqsr_hc_adin, &! Adjoint input : after temperature
530	& zqsr_hc_b_tlout, &! qsr_hc_bngent output: after temperature
531	& zqsr_hc_b_adin, &! Adjoint input : after temperature
532	& zetot3_tlin, &! Tangent input
533	& zetot3_adout, &! Adjoint output
534	& zta, & ! temporary after temperature
535	& zqsr_hc, & ! temporary after temperature
536	& zqsr_hc_b, & ! temporary after temperature
537	& zetot3 ! temporary
538	REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
539	& zqsr_tlin, &! Tangent input : solar radiation (w/m2)
540	& zqsr_adout, &! Adjoint output: solar radiation (w/m2)
541	& zqsr ! temporary solar radiation (w/m2)
542	REAL(KIND=wp) :: &
543	& zsp1, & ! scalar product involving the tangent routine
544	& zsp2, & ! scalar product involving the adjoint routine
545	& zsp2_1, & ! scalar product involving the adjoint routine
546	& zsp2_2, & ! scalar product involving the adjoint routine
547	& zsp2_3 ! scalar product involving the adjoint routine
548	CHARACTER(LEN=14) :: &
549	& cl_name
550
551	ALLOCATE( &
552	& zta_tlin(jpi,jpj,jpk), &
553	& zta_tlout(jpi,jpj,jpk), &
554	& zta_adout(jpi,jpj,jpk), &
555	& zta_adin(jpi,jpj,jpk), &
556	& zta(jpi,jpj,jpk), &
557	& zqsr_hc_tlin(jpi,jpj,jpk), &
558	& zqsr_hc_tlout(jpi,jpj,jpk), &
559	& zqsr_hc_adout(jpi,jpj,jpk), &
560	& zqsr_hc_adin(jpi,jpj,jpk), &
561	& zqsr_hc(jpi,jpj,jpk), &
562	& zqsr_hc_b_tlout(jpi,jpj,jpk), &
563	& zqsr_hc_b_adin(jpi,jpj,jpk), &
564	& zqsr_hc_b(jpi,jpj,jpk), &
565	& zqsr_tlin(jpi,jpj), &
566	& zqsr_adout(jpi,jpj), &
567	& zetot3_tlin(jpi,jpj,jpk), &
568	& zetot3_adout(jpi,jpj,jpk), &
569	& zqsr(jpi,jpj), &
570	& zetot3(jpi,jpj,jpk) &
571	& )
572	! Initialize the reference state
573	qsr(:,:) = 1.0_wp ! ???
574	!Initialize etot3 to non-zero value until traj(nit000-1) is fixed
575	etot3(:,:,1) = 2.e-8 ; etot3(:,:,2) = 1.5e-9; etot3(:,:,3) = 8.5e-10
576	etot3(:,:,4) = 5.4e-10 ; etot3(:,:,5) = 3.5e-10; etot3(:,:,6:jpk) = 0.0_wp
577	! Initialize random field standard deviations
578	!=============================================================
579	! 1) dx = ( T ) and dy = ( T )
580	!=============================================================
581
582	!--------------------------------------------------------------------
583	! Reset the tangent and adjoint variables
584	!--------------------------------------------------------------------
585	zta_tlin(:,:,:) = 0.0_wp
586	zta_tlout(:,:,:) = 0.0_wp
587	zta_adout(:,:,:) = 0.0_wp
588	zta_adin(:,:,:) = 0.0_wp
589	zqsr_hc_tlin(:,:,:) = 0.0_wp
590	zqsr_hc_tlout(:,:,:) = 0.0_wp
591	zqsr_hc_adout(:,:,:) = 0.0_wp
592	zqsr_hc_adin(:,:,:) = 0.0_wp
593	zqsr_hc_b_tlout(:,:,:) = 0.0_wp
594	zqsr_hc_b_adin(:,:,:) = 0.0_wp
595	zqsr_adout(:,:) = 0.0_wp
596	zqsr_tlin(:,:) = 0.0_wp
597	zetot3_tlin(:,:,:) = 0.0_wp
598	zetot3_adout(:,:,:) = 0.0_wp
599	tsa_ad(:,:,:,jp_tem) = 0.0_wp
600	qsr_ad(:,:) = 0.0_wp
601	qsr_hc_ad(:,:,:) = 0.0_wp
602	qsr_hc_b_ad(:,:,:) = 0.0_wp
603	etot3_ad(:,:,:) = 0.0_wp
604
605	CALL grid_random( zqsr , 'T', 0.0_wp, stdqsr )
606	CALL grid_random( zqsr_hc, 'T', 0.0_wp, stdqsr )
607	CALL grid_random( zta , 'T', 0.0_wp, stdt )
608	CALL grid_random( zetot3 , 'T', 0.0_wp, stdt )
609	DO jk = 1, jpk
610	DO jj = nldj, nlej
611	DO ji = nldi, nlei
612	zta_tlin(ji,jj,jk) = zta(ji,jj,jk)
613	END DO
614	END DO
615	END DO
616	DO jk = 1, jpk
617	DO jj = nldj, nlej
618	DO ji = nldi, nlei
619	zqsr_hc_tlin(ji,jj,jk) = zqsr_hc(ji,jj,jk)
620	END DO
621	END DO
622	END DO
623	DO jk = 1, jpk
624	DO jj = nldj, nlej
625	DO ji = nldi, nlei
626	zetot3_tlin(ji,jj,jk) = zetot3(ji,jj,jk)
627	END DO
628	END DO
629	END DO
630	DO jj = nldj, nlej
631	DO ji = nldi, nlei
632	zqsr_tlin(ji,jj) = zqsr(ji,jj)
633	END DO
634	END DO
635	! Test for time steps nit000 and nit000 + 1 (the matrix changes)
636	DO jstp = nit000, nit000 + 1
637	!--------------------------------------------------------------------
638	! Call the tangent routine: dy = L dx
639	!--------------------------------------------------------------------
640
641	tsa_tl(:,:,:,jp_tem) = zta_tlin(:,:,:)
642	etot3_tl(:,:,:) = zetot3_tlin(:,:,:)
643	qsr_tl(:,:) = zqsr_tlin(:,:)
644	qsr_hc_tl(:,:,:) = zqsr_hc_tlin(:,:,:)
645
646	CALL tra_qsr_tan( jstp )
647
648	zta_tlout(:,:,:) = tsa_tl(:,:,:,jp_tem)
649	zqsr_hc_tlout(:,:,:) = qsr_hc_tl(:,:,:)
650	zqsr_hc_b_tlout(:,:,:) = qsr_hc_b_tl(:,:,:)
651
652	!--------------------------------------------------------------------
653	! Initialize the adjoint variables: dy^* = W dy
654	!--------------------------------------------------------------------
655
656	DO jk = 1, jpk
657	DO jj = nldj, nlej
658	DO ji = nldi, nlei
659	zta_adin(ji,jj,jk) = zta_tlout(ji,jj,jk) &
660	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
661	& * tmask(ji,jj,jk)
662	END DO
663	END DO
664	END DO
665	DO jk = 1, jpk
666	DO jj = nldj, nlej
667	DO ji = nldi, nlei
668	zqsr_hc_adin(ji,jj,jk) = zqsr_hc_tlout(ji,jj,jk) &
669	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
670	& * tmask(ji,jj,jk)
671	END DO
672	END DO
673	END DO
674	DO jk = 1, jpk
675	DO jj = nldj, nlej
676	DO ji = nldi, nlei
677	zqsr_hc_b_adin(ji,jj,jk) = zqsr_hc_b_tlout(ji,jj,jk) &
678	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
679	& * tmask(ji,jj,jk)
680	END DO
681	END DO
682	END DO
683
684	!--------------------------------------------------------------------
685	! Compute the scalar product: ( L dx )^T W dy
686	!--------------------------------------------------------------------
687
688	zsp1 = DOT_PRODUCT( zta_tlout, zta_adin ) &
689	& + DOT_PRODUCT( zqsr_hc_tlout, zqsr_hc_adin ) &
690	& + DOT_PRODUCT( zqsr_hc_b_tlout, zqsr_hc_b_adin )
691
692	!--------------------------------------------------------------------
693	! Call the adjoint routine: dx^* = L^T dy^*
694	!--------------------------------------------------------------------
695
696	etot3_ad(:,:,:) = 0.0_wp
697	qsr_ad(:,:) = 0.0_wp
698	tsa_ad(:,:,:,jp_tem) = zta_adin(:,:,:)
699	qsr_hc_ad(:,:,:) = zqsr_hc_adin(:,:,:)
700	qsr_hc_b_ad(:,:,:) = zqsr_hc_b_adin(:,:,:)
701
702	CALL tra_qsr_adj( jstp )
703
704	zta_adout(:,:,:) = tsa_ad(:,:,:,jp_tem)
705	zetot3_adout(:,:,:) = etot3_ad(:,:,:)
706	zqsr_adout(:,:) = qsr_ad(:,:)
707	zqsr_hc_adout(:,:,:) = qsr_hc_ad(:,:,:)
708
709	!--------------------------------------------------------------------
710	! Compute the scalar product: dx^T L^T W dy
711	!--------------------------------------------------------------------
712
713	zsp2_1 = DOT_PRODUCT( zta_tlin , zta_adout )
714	zsp2_1 = zsp2_1 + DOT_PRODUCT( zqsr_hc_tlin , zqsr_hc_adout )
715	zsp2_2 = DOT_PRODUCT( zqsr_tlin , zqsr_adout )
716	zsp2_3 = DOT_PRODUCT( zetot3_tlin , zetot3_adout )
717
718	zsp2 = zsp2_1 + zsp2_2 + zsp2_3
719
720	! Compare the scalar products
721
722	! 14 char: '12345678901234'
723	IF (jstp == nit000) THEN
724	cl_name = 'tra_qsr_adj 1'
725	ELSE
726	cl_name = 'tra_qsr_adj 2'
727	END IF
728	CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
729	END DO
730
731	DEALLOCATE( &
732	& zta_tlin, &
733	& zta_tlout, &
734	& zta_adout, &
735	& zta_adin, &
736	& zta, &
737	& zqsr_hc_tlin, &
738	& zqsr_hc_tlout, &
739	& zqsr_hc_adout, &
740	& zqsr_hc_adin, &
741	& zqsr_hc, &
742	& zqsr_hc_b_tlout, &
743	& zqsr_hc_b_adin, &
744	& zqsr_hc_b, &
745	& zqsr_adout, &
746	& zqsr_tlin, &
747	& zqsr &
748	& )
749
750	!
751	END SUBROUTINE tra_qsr_adj_tst
752	SUBROUTINE tra_qsr_init_tam
753	!!----------------------------------------------------------------------
754	!! * ROUTINE tra_qsr_init_tan *
755	!!
756	!! ** Purpose : Initialization for the penetrative solar radiation
757	!!
758	!! ** Method : The profile of solar radiation within the ocean is set
759	!! from two length scale of penetration (rn_si0,rn_si1) and a ratio
760	!! (rn_abs). These parameters are read in the namtra_qsr namelist. The
761	!! default values correspond to clear water (type I in Jerlov'
762	!! (1968) classification.
763	!! called by tra_qsr at the first timestep (nit000)
764	!!
765	!! ** Action : - initialize rn_si0, rn_si1 and rn_abs
766	!!
767	!! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
768	!!----------------------------------------------------------------------
769
770	IF( ln_traqsr ) THEN ! Initialisation of Light Penetration !
771	! ! ===================================== !
772	!
773	! ! ---------------------------------- !
774	IF( ln_qsr_rgb ) THEN ! Red-Green-Blue light penetration !
775	! ! ---------------------------------- !
776	etot3_tl(:,:,:) = 0.0_wp
777	etot3_ad(:,:,:) = 0.0_wp
778	!
779	ENDIF
780	! ! ---------------------------------- !
781	IF( ln_qsr_2bd ) THEN ! 2 bands light penetration !
782	! ! ---------------------------------- !
783	etot3_tl(:,:,:) = 0.0_wp
784	etot3_ad(:,:,:) = 0.0_wp
785	!
786	ENDIF
787	!
788	ENDIF
789	!
790	END SUBROUTINE tra_qsr_init_tam
791
792	!!======================================================================
793	#endif
794	END MODULE traqsr_tam

Note: See TracBrowser for help on using the repository browser.

Download in other formats: