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

Last change on this file since 3611 was 3611, checked in by pabouttier, 11 years ago
Add TAM code and ORCA2_TAM configuration - see Ticket #1007
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) not implemented in TAM yet')
213	ELSE !* Constant Chlorophyll
214	DO jk = 1, nksr
215	qsr_hc_tl(:,:,jk) = etot3_tl(:,:,jk) * qsr(:,:) + etot3(:,:,jk) * qsr_tl(:,:)
216	END DO
217	ENDIF
218
219	ENDIF
220	! ! ------------------------- !
221	IF( ln_qsr_2bd ) THEN ! 2 band light penetration !
222	! ! ------------------------- !
223	!
224	IF( lk_vvl ) THEN !* variable volume
225	!zz0 = rn_abs * ro0cpr
226	!zz1 = ( 1. - rn_abs ) * ro0cpr
227	!DO jk = 1, nksr ! solar heat absorbed at T-point in the top 400m
228	!DO jj = 1, jpj
229	!DO ji = 1, jpi
230	!zc0 = zz0 * EXP( -fsdepw(ji,jj,jk )xsi0r ) + zz1 EXP( -fsdepw(ji,jj,jk )*xsi1r )
231	!zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)xsi0r ) + zz1 EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
232	!qsr_hc_tl(ji,jj,jk) = qsr_tl(ji,jj) * ( zc0tmask(ji,jj,jk) - zc1tmask(ji,jj,jk+1) )
233	!END DO
234	!END DO
235	!END DO
236	CALL ctl_stop('tra_qsr_tan: key_vvl or chlorophyll management not implemented in TAM yet')
237	ELSE !* constant volume: coef. computed one for all
238	DO jk = 1, nksr
239	DO jj = 2, jpjm1
240	DO ji = fs_2, fs_jpim1 ! vector opt.
241	qsr_hc_tl(ji,jj,jk) = etot3_tl(ji,jj,jk) * qsr(ji,jj) + etot3(ji,jj,jk) * qsr_tl(ji,jj)
242	END DO
243	END DO
244	END DO
245	!
246	ENDIF
247	!
248	ENDIF
249	DO jk = 1, nksr
250	DO jj = 2, jpjm1
251	DO ji = fs_2, fs_jpim1 ! vector opt.
252	z1_e3t = zfact / fse3t(ji,jj,jk)
253	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
254	END DO
255	END DO
256	END DO
257	!
258	ENDIF
259	!
260	CALL wrk_dealloc( jpi, jpj, zekb, zekg, zekr )
261	CALL wrk_dealloc( jpi, jpj, jpk, ze0tl, ze1tl, ze2tl, ze3tl, zeatl )
262	!
263	IF( nn_timing == 1 ) CALL timing_stop('tra_qsr_tan')
264	!
265	END SUBROUTINE tra_qsr_tan
266	SUBROUTINE tra_qsr_adj( kt )
267	!!----------------------------------------------------------------------
268	!! * ROUTINE tra_qsr_adj *
269	!!
270	!! ** Purpose : Compute the temperature trend due to the solar radiation
271	!! penetration and add it to the general temperature trend.
272	!!
273	!! ** Method : The profile of the solar radiation within the ocean is defined
274	!! through 2 wavebands (rn_si0,rn_si1) or 3 wavebands (RGB) and a ratio rn_abs
275	!! Considering the 2 wavebands case:
276	!! I(k) = Qsr( rn_absEXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) )
277	!! The temperature trend associated with the solar radiation penetration
278	!! is given by : zta = 1/e3t dk[ I ] / (rau0*Cp)
279	!! At the bottom, boudary condition for the radiation is no flux :
280	!! all heat which has not been absorbed in the above levels is put
281	!! in the last ocean level.
282	!! In z-coordinate case, the computation is only done down to the
283	!! level where I(k) < 1.e-15 W/m2. In addition, the coefficients
284	!! used for the computation are calculated one for once as they
285	!! depends on k only.
286	!!
287	!! ** Action : - update ta with the penetrative solar radiation trend
288	!!
289	!! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
290	!! Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516.
291	!!----------------------------------------------------------------------
292	!!
293	INTEGER, INTENT(in) :: kt ! ocean time-step
294	!
295	!!
296	INTEGER :: ji, jj, jk ! dummy loop indexes
297	INTEGER :: irgb ! temporary integers
298	REAL(wp) :: zchl, zcoef, zfact, z1_e3t ! temporary scalars
299	REAL(wp) :: zc0, zc1, zc2, zc3, zz0, zz1 ! - -
300	REAL(wp), POINTER, DIMENSION(:,:) :: zekb, zekg, zekr ! 2D workspace
301	REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0ad, ze1ad , ze2ad, ze3ad, zeaad ! 3D workspace
302	!!----------------------------------------------------------------------
303	!
304	IF( nn_timing == 1 ) CALL timing_start('tra_qsr_adj')
305	!
306	CALL wrk_alloc( jpi, jpj, zekb, zekg, zekr )
307	CALL wrk_alloc( jpi, jpj, jpk, ze0ad, ze1ad, ze2ad, ze3ad, zeaad )
308	!
309	IF( kt == nitend ) THEN
310	IF(lwp) WRITE(numout,*)
311	IF(lwp) WRITE(numout,*) 'tra_qsr_adj : penetration of the surface solar radiation'
312	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
313	IF( .NOT.ln_traqsr ) RETURN
314	ENDIF
315	! Set before qsr tracer content field
316	! ***********************************
317	IF( kt == nit000 ) THEN ! Set the forcing field at nit000 - 1
318	! ! -----------------------------------
319	IF( ln_rstart ) THEN !.AND. & ! Restart: read in restart file
320	!& iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN
321	!IF(lwp) WRITE(numout,*) ' nit000-1 qsr tracer content forcing field red in the restart file'
322	zfact = 0.5e0
323	!CALL iom_get( numror, jpdom_autoglo, 'qsr_hc_b', qsr_hc_b_ad ) ! before heat content trend due to Qsr flux
324	ELSE ! No restart or restart not found: Euler forward time stepping
325	zfact = 1.e0
326	ENDIF
327	ELSE ! Swap of forcing field
328	! ! ---------------------
329	zfact = 0.5e0
330	ENDIF
331	! ! ============================================== !
332	IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN ! bio-model fluxes : all vertical coordinates !
333	! ! ============================================== !
334	DO jk = jpkm1, 1, -1
335	DO jj = 2, jpjm1
336	DO ji = fs_2, fs_jpim1 ! vector opt.
337	z1_e3t = zfact / fse3t(ji,jj,jk)
338	qsr_hc_b_ad(ji,jj,jk) = qsr_hc_b_ad(ji,jj,jk) + tsa_ad(ji,jj,jk,jp_tem) * z1_e3t
339	qsr_hc_ad(ji,jj,jk) = qsr_hc_ad(ji,jj,jk) + tsa_ad(ji,jj,jk,jp_tem) * z1_e3t
340	END DO
341	END DO
342	END DO
343	DO jk = jpkm1, 1, -1
344	etot3_ad(:,:,jk) = etot3_ad(:,:,jk) + ro0cpr * qsr_hc_ad(:,:,jk)
345	etot3_ad(:,:,jk+1) = etot3_ad(:,:,jk+1) - ro0cpr * qsr_hc_ad(:,:,jk)
346	END DO
347	! ! ============================================== !
348	ELSE ! Ocean alone :
349	! ! ============================================== !
350	!
351	DO jk = 1, nksr
352	DO jj = 2, jpjm1
353	DO ji = fs_2, fs_jpim1 ! vector opt.
354	z1_e3t = zfact / fse3t(ji,jj,jk)
355	qsr_hc_b_ad(ji,jj,jk) = qsr_hc_b_ad(ji,jj,jk) + tsa_ad(ji,jj,jk,jp_tem) * z1_e3t
356	qsr_hc_ad(ji,jj,jk) = qsr_hc_ad(ji,jj,jk) + tsa_ad(ji,jj,jk,jp_tem) * z1_e3t
357	END DO
358	END DO
359	END DO
360	! ! ------------------------- !
361	IF( ln_qsr_2bd ) THEN ! 2 band light penetration !
362	! ! ------------------------- !
363	!
364	IF( lk_vvl ) THEN !* variable volume
365	!zz0 = rn_abs * ro0cpr
366	!zz1 = ( 1. - rn_abs ) * ro0cpr
367	!DO jk = nksr, 1, -1 ! solar heat absorbed at T-point in the top 400m
368	!DO jj = 1, jpj
369	!DO ji = 1, jpi
370	!zc0 = zz0 * EXP( -fsdepw(ji,jj,jk )xsi0r ) + zz1 EXP( -fsdepw(ji,jj,jk )*xsi1r )
371	!zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)xsi0r ) + zz1 EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
372	!qsr_ad(ji,jj) = qsr_hc_ad(ji,jj) * ( zc0tmask(ji,jj,jk) - zc1tmask(ji,jj,jk+1) )
373	!END DO
374	!END DO
375	!END DO
376	CALL ctl_stop('tra_qsr_adj: key_vvl or chlorophyll management not implemented in TAM yet')
377	ELSE !* constant volume: coef. computed one for all
378	DO jk = 1, nksr
379	DO jj = 2, jpjm1
380	DO ji = fs_2, fs_jpim1 ! vector opt.
381	etot3_ad(ji,jj,jk) = etot3_ad(ji,jj,jk) + qsr(ji,jj) * qsr_hc_ad(ji,jj,jk)
382	qsr_ad(ji,jj) = qsr_ad(ji,jj) + etot3(ji,jj,jk) * qsr_hc_ad(ji,jj,jk)
383	qsr_hc_ad(ji,jj,jk) = 0._wp
384	END DO
385	END DO
386	END DO
387	!
388	ENDIF
389	!
390	ENDIF
391	!
392	! ! ------------------------- !
393	IF( ln_qsr_rgb) THEN ! R-G-B light penetration !
394	! ! ------------------------- !
395	! Set chlorophyl concentration
396	IF( nn_chldta == 1 .OR. lk_vvl ) THEN !* Variable Chlorophyll or ocean volume
397	!!!
398	!!IF( nn_chldta ==1 ) THEN !* Variable Chlorophyll
399	!!zc0ad = 0.0_wp; zc1ad = 0.0_wp; zc2ad = 0.0_wp; zc3ad = 0.0_wp
400	!!ze0ad(:,:,:) = 0.0_wp; ze1ad(:,:,:) = 0.0_wp; ze2ad(:,:,:) = 0.0_wp; ze3ad(:,:,:) = 0.0_wp
401	!!zeaad(:,:,:) = 0.0_wp
402	!!!
403	!!CALL fld_read( kt, 1, sf_chl ) ! Read Chl data and provides it at the current time step
404	!!!
405	!!!CDIR COLLAPSE
406	!!!CDIR NOVERRCHK
407	!!DO jj = 1, jpj ! Separation in R-G-B depending of the surface Chl
408	!!!CDIR NOVERRCHK
409	!!DO ji = 1, jpi
410	!!zchl = MIN( 10.0_wp , MAX( 0.03_wp, sf_chl(1)%fnow(ji,jj) ) )
411	!!irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
412	!!zekb(ji,jj) = rkrgb(1,irgb)
413	!!zekg(ji,jj) = rkrgb(2,irgb)
414	!!zekr(ji,jj) = rkrgb(3,irgb)
415	!!END DO
416	!!END DO
417	!ELSE
418	!zchl = 0.05 ! constant chlorophyll
419	!irgb = NINT( 41 + 20.*LOG10( zchl ) + 1.e-15 )
420	!zekb(:,:) = rkrgb(1,irgb) ! Separation in R-G-B depending of the chlorophyll
421	!zekg(:,:) = rkrgb(2,irgb)
422	!zekr(:,:) = rkrgb(3,irgb)
423	!ENDIF
424	!!
425	!zcoef = ( 1.0_wp - rn_abs ) / 3.0_wp
426
427	!zeaad(:,:,nksr+1:jpk) = 0.0_wp ! below 400m set to zero
428	!!
429	!DO jk = 1, nksr ! compute and add qsr trend to ta
430	!zeaad(:,:,jk ) = ro0cpr * qsr_hc_ad(:,:,jk)
431	!zeaad(:,:,jk+1) = - ro0cpr * qsr_hc_ad(:,:,jk)
432	!END DO
433	!!
434	!DO jk = nksr+1, 2, -1
435	!!CDIR NOVERRCHK
436	!DO jj = 1, jpj
437	!!CDIR NOVERRCHK
438	!DO ji = 1, jpi
439	!zc0ad = zc0ad + zeaad(ji,jj,jk) * tmask(ji,jj,jk)
440	!zc1ad = zc1ad + zeaad(ji,jj,jk) * tmask(ji,jj,jk)
441	!zc2ad = zc2ad + zeaad(ji,jj,jk) * tmask(ji,jj,jk)
442	!zc3ad = zc3ad + zeaad(ji,jj,jk) * tmask(ji,jj,jk)
443	!zeaad(ji,jj,jk) = 0.0_wp
444	!zc0ad = zc0ad + ze0ad(ji,jj,jk)
445	!zc1ad = zc1ad + ze1ad(ji,jj,jk)
446	!zc2ad = zc2ad + ze2ad(ji,jj,jk)
447	!zc3ad = zc3ad + ze3ad(ji,jj,jk)
448	!ze0ad(ji,jj,jk) = 0.0_wp
449	!ze1ad(ji,jj,jk) = 0.0_wp
450	!ze2ad(ji,jj,jk) = 0.0_wp
451	!ze3ad(ji,jj,jk) = 0.0_wp
452	!ze0ad(ji,jj,jk-1) = ze0ad(ji,jj,jk-1) + zc0ad * EXP( - fse3t(ji,jj,jk-1) * xsi0r )
453	!ze1ad(ji,jj,jk-1) = ze1ad(ji,jj,jk-1) + zc1ad * EXP( - fse3t(ji,jj,jk-1) * zekb(ji,jj) )
454	!ze2ad(ji,jj,jk-1) = ze2ad(ji,jj,jk-1) + zc2ad * EXP( - fse3t(ji,jj,jk-1) * zekg(ji,jj) )
455	!ze3ad(ji,jj,jk-1) = ze3ad(ji,jj,jk-1) + zc3ad * EXP( - fse3t(ji,jj,jk-1) * zekr(ji,jj) )
456	!zc0ad = 0.0_wp
457	!zc1ad = 0.0_wp
458	!zc2ad = 0.0_wp
459	!zc3ad = 0.0_wp
460	!END DO
461	!END DO
462	!END DO
463	!!
464	!qsr_ad(:,:) = qsr_ad(:,:) + zeaad(:,:,1)
465	!qsr_ad(:,:) = qsr_ad(:,:) + zcoef * ze3ad(:,:,1)
466	!qsr_ad(:,:) = qsr_ad(:,:) + zcoef * ze2ad(:,:,1)
467	!qsr_ad(:,:) = qsr_ad(:,:) + zcoef * ze1ad(:,:,1)
468	!qsr_ad(:,:) = qsr_ad(:,:) + rn_abs * ze0ad(:,:,1)
469	!!
470	CALL ctl_stop('tra_qsr_adj: key_vvl or chlorophyll management not implemented in TAM yet')
471	ELSE !* Constant Chlorophyll
472	DO jk = 1, nksr
473	etot3_ad(:,:,jk) = etot3_ad(:,:,jk) + qsr_hc_ad(:,:,jk) * qsr(:,:)
474	qsr_ad( :,: ) = qsr_ad(:,:) + qsr_hc_ad(:,:,jk) * etot3(:,:,jk)
475	qsr_hc_ad(:,:,jk) = 0._wp
476	END DO
477	ENDIF
478	ENDIF
479	ENDIF
480	IF ( kt /= nit000 ) THEN
481	qsr_hc_ad(:,:,:) = qsr_hc_ad(:,:,:) + qsr_hc_b_ad(:,:,:)
482	ENDIF
483	qsr_hc_b_ad(:,:,:) = 0._wp
484
485	IF( nn_timing == 1 ) CALL timing_stop('tra_qsr_adj')
486
487	END SUBROUTINE tra_qsr_adj
488	SUBROUTINE tra_qsr_adj_tst ( kumadt )
489	!!-----------------------------------------------------------------------
490	!!
491	!! * ROUTINE tra_sbc_adj_tst : TEST OF tra_sbc_adj *
492	!!
493	!! ** Purpose : Test the adjoint routine.
494	!!
495	!! ** Method : Verify the scalar product
496	!!
497	!! ( L dx )^T W dy = dx^T L^T W dy
498	!!
499	!! where L = tangent routine
500	!! L^T = adjoint routine
501	!! W = diagonal matrix of scale factors
502	!! dx = input perturbation (random field)
503	!! dy = L dx
504	!!
505	!! History :
506	!! ! 08-08 (A. Vidard)
507	!!-----------------------------------------------------------------------
508	!! * Modules used
509
510	!! * Arguments
511	INTEGER, INTENT(IN) :: &
512	& kumadt ! Output unit
513
514	INTEGER :: &
515	& jstp, &
516	& ji, & ! dummy loop indices
517	& jj, &
518	& jk
519	!! * Local declarations
520	REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
521	& zta_tlin, &! Tangent input : after temperature
522	& zta_tlout, &! Tangent output: after temperature
523	& zta_adout, &! Adjoint output: after temperature
524	& zta_adin, &! Adjoint input : after temperature
525	& zqsr_hc_tlin, &! qsr_hcngent input : after temperature
526	& zqsr_hc_tlout, &! qsr_hcngent output: after temperature
527	& zqsr_hc_adout, &! Adjoint output: after temperature
528	& zqsr_hc_adin, &! Adjoint input : after temperature
529	& zqsr_hc_b_tlout, &! qsr_hc_bngent output: after temperature
530	& zqsr_hc_b_adin, &! Adjoint input : after temperature
531	& zetot3_tlin, &! Tangent input
532	& zetot3_adout, &! Adjoint output
533	& zta, & ! temporary after temperature
534	& zqsr_hc, & ! temporary after temperature
535	& zqsr_hc_b, & ! temporary after temperature
536	& zetot3 ! temporary
537	REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
538	& zqsr_tlin, &! Tangent input : solar radiation (w/m2)
539	& zqsr_adout, &! Adjoint output: solar radiation (w/m2)
540	& zqsr ! temporary solar radiation (w/m2)
541	REAL(KIND=wp) :: &
542	& zsp1, & ! scalar product involving the tangent routine
543	& zsp2, & ! scalar product involving the adjoint routine
544	& zsp2_1, & ! scalar product involving the adjoint routine
545	& zsp2_2, & ! scalar product involving the adjoint routine
546	& zsp2_3 ! scalar product involving the adjoint routine
547	CHARACTER(LEN=14) :: &
548	& cl_name
549
550	ALLOCATE( &
551	& zta_tlin(jpi,jpj,jpk), &
552	& zta_tlout(jpi,jpj,jpk), &
553	& zta_adout(jpi,jpj,jpk), &
554	& zta_adin(jpi,jpj,jpk), &
555	& zta(jpi,jpj,jpk), &
556	& zqsr_hc_tlin(jpi,jpj,jpk), &
557	& zqsr_hc_tlout(jpi,jpj,jpk), &
558	& zqsr_hc_adout(jpi,jpj,jpk), &
559	& zqsr_hc_adin(jpi,jpj,jpk), &
560	& zqsr_hc(jpi,jpj,jpk), &
561	& zqsr_hc_b_tlout(jpi,jpj,jpk), &
562	& zqsr_hc_b_adin(jpi,jpj,jpk), &
563	& zqsr_hc_b(jpi,jpj,jpk), &
564	& zqsr_tlin(jpi,jpj), &
565	& zqsr_adout(jpi,jpj), &
566	& zetot3_tlin(jpi,jpj,jpk), &
567	& zetot3_adout(jpi,jpj,jpk), &
568	& zqsr(jpi,jpj), &
569	& zetot3(jpi,jpj,jpk) &
570	& )
571	! Initialize the reference state
572	qsr(:,:) = 1.0_wp ! ???
573	!Initialize etot3 to non-zero value until traj(nit000-1) is fixed
574	etot3(:,:,1) = 2.e-8 ; etot3(:,:,2) = 1.5e-9; etot3(:,:,3) = 8.5e-10
575	etot3(:,:,4) = 5.4e-10 ; etot3(:,:,5) = 3.5e-10; etot3(:,:,6:jpk) = 0.0_wp
576	! Initialize random field standard deviations
577	!=============================================================
578	! 1) dx = ( T ) and dy = ( T )
579	!=============================================================
580
581	!--------------------------------------------------------------------
582	! Reset the tangent and adjoint variables
583	!--------------------------------------------------------------------
584	zta_tlin(:,:,:) = 0.0_wp
585	zta_tlout(:,:,:) = 0.0_wp
586	zta_adout(:,:,:) = 0.0_wp
587	zta_adin(:,:,:) = 0.0_wp
588	zqsr_hc_tlin(:,:,:) = 0.0_wp
589	zqsr_hc_tlout(:,:,:) = 0.0_wp
590	zqsr_hc_adout(:,:,:) = 0.0_wp
591	zqsr_hc_adin(:,:,:) = 0.0_wp
592	zqsr_hc_b_tlout(:,:,:) = 0.0_wp
593	zqsr_hc_b_adin(:,:,:) = 0.0_wp
594	zqsr_adout(:,:) = 0.0_wp
595	zqsr_tlin(:,:) = 0.0_wp
596	zetot3_tlin(:,:,:) = 0.0_wp
597	zetot3_adout(:,:,:) = 0.0_wp
598	tsa_ad(:,:,:,jp_tem) = 0.0_wp
599	qsr_ad(:,:) = 0.0_wp
600	qsr_hc_ad(:,:,:) = 0.0_wp
601	qsr_hc_b_ad(:,:,:) = 0.0_wp
602	etot3_ad(:,:,:) = 0.0_wp
603
604	CALL grid_random( zqsr , 'T', 0.0_wp, stdqsr )
605	CALL grid_random( zqsr_hc, 'T', 0.0_wp, stdqsr )
606	CALL grid_random( zta , 'T', 0.0_wp, stdt )
607	CALL grid_random( zetot3 , 'T', 0.0_wp, stdt )
608	DO jk = 1, jpk
609	DO jj = nldj, nlej
610	DO ji = nldi, nlei
611	zta_tlin(ji,jj,jk) = zta(ji,jj,jk)
612	END DO
613	END DO
614	END DO
615	DO jk = 1, jpk
616	DO jj = nldj, nlej
617	DO ji = nldi, nlei
618	zqsr_hc_tlin(ji,jj,jk) = zqsr_hc(ji,jj,jk)
619	END DO
620	END DO
621	END DO
622	DO jk = 1, jpk
623	DO jj = nldj, nlej
624	DO ji = nldi, nlei
625	zetot3_tlin(ji,jj,jk) = zetot3(ji,jj,jk)
626	END DO
627	END DO
628	END DO
629	DO jj = nldj, nlej
630	DO ji = nldi, nlei
631	zqsr_tlin(ji,jj) = zqsr(ji,jj)
632	END DO
633	END DO
634	! Test for time steps nit000 and nit000 + 1 (the matrix changes)
635	DO jstp = nit000, nit000 + 1
636	!--------------------------------------------------------------------
637	! Call the tangent routine: dy = L dx
638	!--------------------------------------------------------------------
639
640	tsa_tl(:,:,:,jp_tem) = zta_tlin(:,:,:)
641	etot3_tl(:,:,:) = zetot3_tlin(:,:,:)
642	qsr_tl(:,:) = zqsr_tlin(:,:)
643	qsr_hc_tl(:,:,:) = zqsr_hc_tlin(:,:,:)
644
645	CALL tra_qsr_tan( jstp )
646
647	zta_tlout(:,:,:) = tsa_tl(:,:,:,jp_tem)
648	zqsr_hc_tlout(:,:,:) = qsr_hc_tl(:,:,:)
649	zqsr_hc_b_tlout(:,:,:) = qsr_hc_b_tl(:,:,:)
650
651	!--------------------------------------------------------------------
652	! Initialize the adjoint variables: dy^* = W dy
653	!--------------------------------------------------------------------
654
655	DO jk = 1, jpk
656	DO jj = nldj, nlej
657	DO ji = nldi, nlei
658	zta_adin(ji,jj,jk) = zta_tlout(ji,jj,jk) &
659	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
660	& * tmask(ji,jj,jk)
661	END DO
662	END DO
663	END DO
664	DO jk = 1, jpk
665	DO jj = nldj, nlej
666	DO ji = nldi, nlei
667	zqsr_hc_adin(ji,jj,jk) = zqsr_hc_tlout(ji,jj,jk) &
668	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
669	& * tmask(ji,jj,jk)
670	END DO
671	END DO
672	END DO
673	DO jk = 1, jpk
674	DO jj = nldj, nlej
675	DO ji = nldi, nlei
676	zqsr_hc_b_adin(ji,jj,jk) = zqsr_hc_b_tlout(ji,jj,jk) &
677	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
678	& * tmask(ji,jj,jk)
679	END DO
680	END DO
681	END DO
682
683	!--------------------------------------------------------------------
684	! Compute the scalar product: ( L dx )^T W dy
685	!--------------------------------------------------------------------
686
687	zsp1 = DOT_PRODUCT( zta_tlout, zta_adin ) &
688	& + DOT_PRODUCT( zqsr_hc_tlout, zqsr_hc_adin ) &
689	& + DOT_PRODUCT( zqsr_hc_b_tlout, zqsr_hc_b_adin )
690
691	!--------------------------------------------------------------------
692	! Call the adjoint routine: dx^* = L^T dy^*
693	!--------------------------------------------------------------------
694
695	etot3_ad(:,:,:) = 0.0_wp
696	qsr_ad(:,:) = 0.0_wp
697	tsa_ad(:,:,:,jp_tem) = zta_adin(:,:,:)
698	qsr_hc_ad(:,:,:) = zqsr_hc_adin(:,:,:)
699	qsr_hc_b_ad(:,:,:) = zqsr_hc_b_adin(:,:,:)
700
701	CALL tra_qsr_adj( jstp )
702
703	zta_adout(:,:,:) = tsa_ad(:,:,:,jp_tem)
704	zetot3_adout(:,:,:) = etot3_ad(:,:,:)
705	zqsr_adout(:,:) = qsr_ad(:,:)
706	zqsr_hc_adout(:,:,:) = qsr_hc_ad(:,:,:)
707
708	!--------------------------------------------------------------------
709	! Compute the scalar product: dx^T L^T W dy
710	!--------------------------------------------------------------------
711
712	zsp2_1 = DOT_PRODUCT( zta_tlin , zta_adout )
713	zsp2_1 = zsp2_1 + DOT_PRODUCT( zqsr_hc_tlin , zqsr_hc_adout )
714	zsp2_2 = DOT_PRODUCT( zqsr_tlin , zqsr_adout )
715	zsp2_3 = DOT_PRODUCT( zetot3_tlin , zetot3_adout )
716
717	zsp2 = zsp2_1 + zsp2_2 + zsp2_3
718
719	! Compare the scalar products
720
721	! 14 char: '12345678901234'
722	IF (jstp == nit000) THEN
723	cl_name = 'tra_qsr_adj 1'
724	ELSE
725	cl_name = 'tra_qsr_adj 2'
726	END IF
727	CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
728	END DO
729
730	DEALLOCATE( &
731	& zta_tlin, &
732	& zta_tlout, &
733	& zta_adout, &
734	& zta_adin, &
735	& zta, &
736	& zqsr_hc_tlin, &
737	& zqsr_hc_tlout, &
738	& zqsr_hc_adout, &
739	& zqsr_hc_adin, &
740	& zqsr_hc, &
741	& zqsr_hc_b_tlout, &
742	& zqsr_hc_b_adin, &
743	& zqsr_hc_b, &
744	& zqsr_adout, &
745	& zqsr_tlin, &
746	& zqsr &
747	& )
748
749	!
750	END SUBROUTINE tra_qsr_adj_tst
751	SUBROUTINE tra_qsr_init_tam
752	!!----------------------------------------------------------------------
753	!! * ROUTINE tra_qsr_init_tan *
754	!!
755	!! ** Purpose : Initialization for the penetrative solar radiation
756	!!
757	!! ** Method : The profile of solar radiation within the ocean is set
758	!! from two length scale of penetration (rn_si0,rn_si1) and a ratio
759	!! (rn_abs). These parameters are read in the namtra_qsr namelist. The
760	!! default values correspond to clear water (type I in Jerlov'
761	!! (1968) classification.
762	!! called by tra_qsr at the first timestep (nit000)
763	!!
764	!! ** Action : - initialize rn_si0, rn_si1 and rn_abs
765	!!
766	!! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
767	!!----------------------------------------------------------------------
768
769	IF( ln_traqsr ) THEN ! Initialisation of Light Penetration !
770	! ! ===================================== !
771	!
772	! ! ---------------------------------- !
773	IF( ln_qsr_rgb ) THEN ! Red-Green-Blue light penetration !
774	! ! ---------------------------------- !
775	etot3_tl(:,:,:) = 0.0_wp
776	etot3_ad(:,:,:) = 0.0_wp
777	!
778	ENDIF
779	! ! ---------------------------------- !
780	IF( ln_qsr_2bd ) THEN ! 2 bands light penetration !
781	! ! ---------------------------------- !
782	etot3_tl(:,:,:) = 0.0_wp
783	etot3_ad(:,:,:) = 0.0_wp
784	!
785	ENDIF
786	!
787	ENDIF
788	!
789	END SUBROUTINE tra_qsr_init_tam
790
791	!!======================================================================
792	#endif
793	END MODULE traqsr_tam

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

Download in other formats: