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

Last change on this file since 3640 was 3640, checked in by pabouttier, 11 years ago
Missing allocation/deallocation in TAM routines - See ticket #1013
Property svn:executable set to ``*
File size: 37.1 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	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
797	END MODULE traqsr_tam

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

Download in other formats: