New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

traqsr.F90 in branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90 @ 4333

Last change on this file since 4333 was 4333, checked in by clem, 10 years ago
remove remaining bugs in LIM3, so that it can run in both regional and global config
Property svn:keywords set to `Id`
File size: 29.9 KB

Rev	Line
[3]	1	MODULE traqsr
	2	!!======================================================================
	3	!! * MODULE traqsr *
	4	!! Ocean physics: solar radiation penetration in the top ocean levels
	5	!!======================================================================
[1423]	6	!! History : OPA ! 1990-10 (B. Blanke) Original code
	7	!! 7.0 ! 1991-11 (G. Madec)
	8	!! ! 1996-01 (G. Madec) s-coordinates
	9	!! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module
	10	!! - ! 2005-11 (G. Madec) zco, zps, sco coordinate
	11	!! 3.2 ! 2009-04 (G. Madec & NEMO team)
[4161]	12	!! 4.0 ! 2012-05 (C. Rousset) store attenuation coef for use in ice model
[3]	13	!!----------------------------------------------------------------------
[503]	14
	15	!!----------------------------------------------------------------------
[3]	16	!! tra_qsr : trend due to the solar radiation penetration
	17	!! tra_qsr_init : solar radiation penetration initialization
	18	!!----------------------------------------------------------------------
	19	USE oce ! ocean dynamics and active tracers
	20	USE dom_oce ! ocean space and time domain
[888]	21	USE sbc_oce ! surface boundary condition: ocean
	22	USE trc_oce ! share SMS/Ocean variables
	23	USE trdmod_oce ! ocean variables trends
[2528]	24	USE trdtra ! ocean active tracers trends
[3]	25	USE in_out_manager ! I/O manager
	26	USE phycst ! physical constants
[258]	27	USE prtctl ! Print control
[1423]	28	USE iom ! I/O manager
	29	USE fldread ! read input fields
[4161]	30	USE restart ! ocean restart
[2715]	31	USE lib_mpp ! MPP library
[3294]	32	USE wrk_nemo ! Memory Allocation
	33	USE timing ! Timing
[4207]	34	USE sbc_ice, ONLY : lk_lim3
[3]	35
	36	IMPLICIT NONE
	37	PRIVATE
	38
[2528]	39	PUBLIC tra_qsr ! routine called by step.F90 (ln_traqsr=T)
	40	PUBLIC tra_qsr_init ! routine called by opa.F90
[3]	41
[4147]	42	! !!* Namelist namtra_qsr: penetrative solar radiation
	43	LOGICAL , PUBLIC :: ln_traqsr !: light absorption (qsr) flag
	44	LOGICAL , PUBLIC :: ln_qsr_rgb !: Red-Green-Blue light absorption flag
	45	LOGICAL , PUBLIC :: ln_qsr_2bd !: 2 band light absorption flag
	46	LOGICAL , PUBLIC :: ln_qsr_bio !: bio-model light absorption flag
[4205]	47	LOGICAL , PUBLIC :: ln_qsr_ice !: light penetration for ice-model LIM3 (clem)
[4147]	48	INTEGER , PUBLIC :: nn_chldta !: use Chlorophyll data (=1) or not (=0)
	49	REAL(wp), PUBLIC :: rn_abs !: fraction absorbed in the very near surface (RGB & 2 bands)
	50	REAL(wp), PUBLIC :: rn_si0 !: very near surface depth of extinction (RGB & 2 bands)
	51	REAL(wp), PUBLIC :: rn_si1 !: deepest depth of extinction (water type I) (2 bands)
[187]	52
[1445]	53	! Module variables
[2528]	54	REAL(wp) :: xsi0r !: inverse of rn_si0
	55	REAL(wp) :: xsi1r !: inverse of rn_si1
[1423]	56	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_chl ! structure of input Chl (file informations, fields read)
[2528]	57	INTEGER, PUBLIC :: nksr ! levels below which the light cannot penetrate ( depth larger than 391 m)
[1445]	58	REAL(wp), DIMENSION(3,61) :: rkrgb !: tabulated attenuation coefficients for RGB absorption
[3]	59
	60	!! * Substitutions
	61	# include "domzgr_substitute.h90"
	62	# include "vectopt_loop_substitute.h90"
	63	!!----------------------------------------------------------------------
[2528]	64	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[888]	65	!! $Id$
[2715]	66	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	67	!!----------------------------------------------------------------------
	68	CONTAINS
	69
	70	SUBROUTINE tra_qsr( kt )
	71	!!----------------------------------------------------------------------
	72	!! * ROUTINE tra_qsr *
	73	!!
	74	!! ** Purpose : Compute the temperature trend due to the solar radiation
	75	!! penetration and add it to the general temperature trend.
	76	!!
[1423]	77	!! ** Method : The profile of the solar radiation within the ocean is defined
	78	!! through 2 wavebands (rn_si0,rn_si1) or 3 wavebands (RGB) and a ratio rn_abs
	79	!! Considering the 2 wavebands case:
	80	!! I(k) = Qsr( rn_absEXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) )
	81	!! The temperature trend associated with the solar radiation penetration
	82	!! is given by : zta = 1/e3t dk[ I ] / (rau0*Cp)
[3]	83	!! At the bottom, boudary condition for the radiation is no flux :
	84	!! all heat which has not been absorbed in the above levels is put
	85	!! in the last ocean level.
	86	!! In z-coordinate case, the computation is only done down to the
	87	!! level where I(k) < 1.e-15 W/m2. In addition, the coefficients
	88	!! used for the computation are calculated one for once as they
	89	!! depends on k only.
	90	!!
	91	!! ** Action : - update ta with the penetrative solar radiation trend
	92	!! - save the trend in ttrd ('key_trdtra')
[1423]	93	!!
	94	!! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
	95	!! Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516.
[503]	96	!!----------------------------------------------------------------------
[2715]	97	!
[503]	98	INTEGER, INTENT(in) :: kt ! ocean time-step
[2715]	99	!
[1423]	100	INTEGER :: ji, jj, jk ! dummy loop indices
[2715]	101	INTEGER :: irgb ! local integers
	102	REAL(wp) :: zchl, zcoef, zfact ! local scalars
[1423]	103	REAL(wp) :: zc0, zc1, zc2, zc3 ! - -
[4161]	104	REAL(wp) :: zzc0, zzc1, zzc2, zzc3 ! - -
[2715]	105	REAL(wp) :: zz0, zz1, z1_e3t ! - -
[3294]	106	REAL(wp), POINTER, DIMENSION(:,: ) :: zekb, zekg, zekr
	107	REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt
[3]	108	!!----------------------------------------------------------------------
[3294]	109	!
	110	IF( nn_timing == 1 ) CALL timing_start('tra_qsr')
	111	!
	112	CALL wrk_alloc( jpi, jpj, zekb, zekg, zekr )
	113	CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )
	114	!
[3]	115	IF( kt == nit000 ) THEN
[503]	116	IF(lwp) WRITE(numout,*)
	117	IF(lwp) WRITE(numout,*) 'tra_qsr : penetration of the surface solar radiation'
	118	IF(lwp) WRITE(numout,*) '~~~~~~~'
[1423]	119	IF( .NOT.ln_traqsr ) RETURN
[3]	120	ENDIF
	121
[503]	122	IF( l_trdtra ) THEN ! Save ta and sa trends
[3294]	123	CALL wrk_alloc( jpi, jpj, jpk, ztrdt )
	124	ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
[216]	125	ENDIF
	126
[2528]	127	! Set before qsr tracer content field
	128	! ***********************************
	129	IF( kt == nit000 ) THEN ! Set the forcing field at nit000 - 1
	130	! ! -----------------------------------
	131	IF( ln_rstart .AND. & ! Restart: read in restart file
	132	& iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN
	133	IF(lwp) WRITE(numout,*) ' nit000-1 qsr tracer content forcing field red in the restart file'
	134	zfact = 0.5e0
	135	CALL iom_get( numror, jpdom_autoglo, 'qsr_hc_b', qsr_hc_b ) ! before heat content trend due to Qsr flux
	136	ELSE ! No restart or restart not found: Euler forward time stepping
	137	zfact = 1.e0
	138	qsr_hc_b(:,:,:) = 0.e0
	139	ENDIF
	140	ELSE ! Swap of forcing field
	141	! ! ---------------------
	142	zfact = 0.5e0
	143	qsr_hc_b(:,:,:) = qsr_hc(:,:,:)
	144	ENDIF
	145	! Compute now qsr tracer content field
	146	! ************************************
[1423]	147
	148	! ! ============================================== !
[1445]	149	IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN ! bio-model fluxes : all vertical coordinates !
[1423]	150	! ! ============================================== !
	151	DO jk = 1, jpkm1
[3625]	152	qsr_hc(:,:,jk) = r1_rau0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) )
[2528]	153	END DO
	154	! Add to the general trend
	155	DO jk = 1, jpkm1
	156	DO jj = 2, jpjm1
[3]	157	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	158	z1_e3t = zfact / fse3t(ji,jj,jk)
	159	tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) * z1_e3t
[3]	160	END DO
	161	END DO
[1423]	162	END DO
[1756]	163	CALL iom_put( 'qsr3d', etot3 ) ! Shortwave Radiation 3D distribution
[4161]	164	! clem: store attenuation coefficient of the first ocean level
[4206]	165	IF ( lk_lim3 .AND. ln_qsr_ice ) THEN
[4161]	166	DO jj = 1, jpj
	167	DO ji = 1, jpi
	168	IF ( qsr(ji,jj) /= 0._wp ) THEN
	169	oatte(ji,jj) = ( qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) ) )
	170	iatte(ji,jj) = oatte(ji,jj)
	171	ENDIF
	172	END DO
	173	END DO
	174	ENDIF
[1423]	175	! ! ============================================== !
	176	ELSE ! Ocean alone :
	177	! ! ============================================== !
	178	!
	179	! ! ------------------------- !
	180	IF( ln_qsr_rgb) THEN ! R-G-B light penetration !
	181	! ! ------------------------- !
	182	! Set chlorophyl concentration
[2528]	183	IF( nn_chldta == 1 .OR. lk_vvl ) THEN !* Variable Chlorophyll or ocean volume
[1423]	184	!
[2528]	185	IF( nn_chldta == 1 ) THEN !* Variable Chlorophyll
	186	!
	187	CALL fld_read( kt, 1, sf_chl ) ! Read Chl data and provides it at the current time step
	188	!
[1423]	189	!CDIR COLLAPSE
	190	!CDIR NOVERRCHK
[2528]	191	DO jj = 1, jpj ! Separation in R-G-B depending of the surface Chl
[1423]	192	!CDIR NOVERRCHK
[2528]	193	DO ji = 1, jpi
	194	zchl = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) )
	195	irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
	196	zekb(ji,jj) = rkrgb(1,irgb)
	197	zekg(ji,jj) = rkrgb(2,irgb)
	198	zekr(ji,jj) = rkrgb(3,irgb)
	199	END DO
[187]	200	END DO
[2528]	201	ELSE ! Variable ocean volume but constant chrlorophyll
	202	zchl = 0.05 ! constant chlorophyll
	203	irgb = NINT( 41 + 20.*LOG10( zchl ) + 1.e-15 )
	204	zekb(:,:) = rkrgb(1,irgb) ! Separation in R-G-B depending of the chlorophyll
	205	zekg(:,:) = rkrgb(2,irgb)
	206	zekr(:,:) = rkrgb(3,irgb)
	207	ENDIF
[1423]	208	!
	209	zcoef = ( 1. - rn_abs ) / 3.e0 ! equi-partition in R-G-B
	210	ze0(:,:,1) = rn_abs * qsr(:,:)
	211	ze1(:,:,1) = zcoef * qsr(:,:)
	212	ze2(:,:,1) = zcoef * qsr(:,:)
	213	ze3(:,:,1) = zcoef * qsr(:,:)
	214	zea(:,:,1) = qsr(:,:)
	215	!
	216	DO jk = 2, nksr+1
	217	!CDIR NOVERRCHK
	218	DO jj = 1, jpj
	219	!CDIR NOVERRCHK
	220	DO ji = 1, jpi
[2528]	221	zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r )
[1423]	222	zc1 = ze1(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekb(ji,jj) )
	223	zc2 = ze2(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekg(ji,jj) )
	224	zc3 = ze3(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekr(ji,jj) )
	225	ze0(ji,jj,jk) = zc0
	226	ze1(ji,jj,jk) = zc1
	227	ze2(ji,jj,jk) = zc2
	228	ze3(ji,jj,jk) = zc3
	229	zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,jk)
	230	END DO
	231	END DO
	232	END DO
[4161]	233	! clem: store attenuation coefficient of the first ocean level
[4206]	234	IF ( lk_lim3 .AND. ln_qsr_ice ) THEN
[4161]	235	DO jj = 1, jpj
	236	DO ji = 1, jpi
	237	zzc0 = rn_abs * EXP( - fse3t(ji,jj,1) * xsi0r )
	238	zzc1 = zcoef * EXP( - fse3t(ji,jj,1) * zekb(ji,jj) )
	239	zzc2 = zcoef * EXP( - fse3t(ji,jj,1) * zekg(ji,jj) )
	240	zzc3 = zcoef * EXP( - fse3t(ji,jj,1) * zekr(ji,jj) )
	241	oatte(ji,jj) = 1.0 - ( zzc0 + zzc1 + zzc2 + zzc3 ) * tmask(ji,jj,2)
	242	iatte(ji,jj) = 1.0 - ( zzc0 + zzc1 + zcoef + zcoef ) * tmask(ji,jj,2)
	243	END DO
	244	END DO
	245	ENDIF
[1423]	246	!
	247	DO jk = 1, nksr ! compute and add qsr trend to ta
[3625]	248	qsr_hc(:,:,jk) = r1_rau0_rcp * ( zea(:,:,jk) - zea(:,:,jk+1) )
[1423]	249	END DO
[1756]	250	zea(:,:,nksr+1:jpk) = 0.e0 ! below 400m set to zero
	251	CALL iom_put( 'qsr3d', zea ) ! Shortwave Radiation 3D distribution
[1423]	252	!
	253	ELSE !* Constant Chlorophyll
	254	DO jk = 1, nksr
[2528]	255	qsr_hc(:,:,jk) = etot3(:,:,jk) * qsr(:,:)
[1423]	256	END DO
[4161]	257	! clem: store attenuation coefficient of the first ocean level
[4206]	258	IF ( lk_lim3 .AND. ln_qsr_ice ) THEN
[4161]	259	oatte(:,:) = etot3(:,:,1) / r1_rau0_rcp
	260	iatte(:,:) = oatte(:,:)
	261	ENDIF
	262	ENDIF
[1423]	263
[1448]	264	ENDIF
	265	! ! ------------------------- !
	266	IF( ln_qsr_2bd ) THEN ! 2 band light penetration !
[1423]	267	! ! ------------------------- !
	268	!
[2528]	269	IF( lk_vvl ) THEN !* variable volume
[3625]	270	zz0 = rn_abs * r1_rau0_rcp
	271	zz1 = ( 1. - rn_abs ) * r1_rau0_rcp
[2528]	272	DO jk = 1, nksr ! solar heat absorbed at T-point in the top 400m
[3294]	273	DO jj = 1, jpj
	274	DO ji = 1, jpi
[2528]	275	zc0 = zz0 * EXP( -fsdepw(ji,jj,jk )xsi0r ) + zz1 EXP( -fsdepw(ji,jj,jk )*xsi1r )
	276	zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)xsi0r ) + zz1 EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
	277	qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0tmask(ji,jj,jk) - zc1tmask(ji,jj,jk+1) )
	278	END DO
[187]	279	END DO
	280	END DO
[4161]	281	! clem: store attenuation coefficient of the first ocean level
[4206]	282	IF ( lk_lim3 .AND. ln_qsr_ice ) THEN
[4161]	283	DO jj = 1, jpj
	284	DO ji = 1, jpi
	285	zc0 = zz0 * EXP( -fsdepw(ji,jj,1)xsi0r ) + zz1 EXP( -fsdepw(ji,jj,1)*xsi1r )
	286	zc1 = zz0 * EXP( -fsdepw(ji,jj,2)xsi0r ) + zz1 EXP( -fsdepw(ji,jj,2)*xsi1r )
	287	oatte(ji,jj) = ( zc0tmask(ji,jj,1) - zc1tmask(ji,jj,2) ) / r1_rau0_rcp
	288	iatte(ji,jj) = oatte(ji,jj)
	289	END DO
	290	END DO
	291	ENDIF
[2528]	292	ELSE !* constant volume: coef. computed one for all
	293	DO jk = 1, nksr
	294	DO jj = 2, jpjm1
	295	DO ji = fs_2, fs_jpim1 ! vector opt.
	296	qsr_hc(ji,jj,jk) = etot3(ji,jj,jk) * qsr(ji,jj)
	297	END DO
	298	END DO
	299	END DO
[4161]	300	! clem: store attenuation coefficient of the first ocean level
[4206]	301	IF ( lk_lim3 .AND. ln_qsr_ice ) THEN
[4161]	302	oatte(:,:) = etot3(:,:,1) / r1_rau0_rcp
	303	iatte(:,:) = oatte(:,:)
	304	ENDIF
[2528]	305	!
	306	ENDIF
[1423]	307	!
[187]	308	ENDIF
	309	!
[2528]	310	! Add to the general trend
	311	DO jk = 1, nksr
	312	DO jj = 2, jpjm1
	313	DO ji = fs_2, fs_jpim1 ! vector opt.
	314	z1_e3t = zfact / fse3t(ji,jj,jk)
	315	tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) * z1_e3t
	316	END DO
	317	END DO
	318	END DO
	319	!
[3]	320	ENDIF
[2528]	321	!
	322	IF( lrst_oce ) THEN ! Write in the ocean restart file
	323	! *******************************
	324	IF(lwp) WRITE(numout,*)
	325	IF(lwp) WRITE(numout,*) 'qsr tracer content forcing field written in ocean restart file ', &
	326	& 'at it= ', kt,' date= ', ndastp
	327	IF(lwp) WRITE(numout,*) '~~~~'
	328	CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b', qsr_hc )
	329	!
	330	ENDIF
[3]	331
[503]	332	IF( l_trdtra ) THEN ! qsr tracers trends saved for diagnostics
[2528]	333	ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
	334	CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_qsr, ztrdt )
[3294]	335	CALL wrk_dealloc( jpi, jpj, jpk, ztrdt )
[3]	336	ENDIF
[457]	337	! ! print mean trends (used for debugging)
[2528]	338	IF(ln_ctl) CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' qsr - Ta: ', mask1=tmask, clinfo3='tra-ta' )
[503]	339	!
[3294]	340	CALL wrk_dealloc( jpi, jpj, zekb, zekg, zekr )
	341	CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )
[2715]	342	!
[3294]	343	IF( nn_timing == 1 ) CALL timing_stop('tra_qsr')
	344	!
[3]	345	END SUBROUTINE tra_qsr
	346
	347
	348	SUBROUTINE tra_qsr_init
	349	!!----------------------------------------------------------------------
	350	!! * ROUTINE tra_qsr_init *
	351	!!
	352	!! ** Purpose : Initialization for the penetrative solar radiation
	353	!!
	354	!! ** Method : The profile of solar radiation within the ocean is set
[1423]	355	!! from two length scale of penetration (rn_si0,rn_si1) and a ratio
[1601]	356	!! (rn_abs). These parameters are read in the namtra_qsr namelist. The
[3]	357	!! default values correspond to clear water (type I in Jerlov'
	358	!! (1968) classification.
	359	!! called by tra_qsr at the first timestep (nit000)
	360	!!
[1423]	361	!! ** Action : - initialize rn_si0, rn_si1 and rn_abs
[3]	362	!!
[503]	363	!! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
[3]	364	!!----------------------------------------------------------------------
[2715]	365	!
[4147]	366	INTEGER :: ji, jj, jk ! dummy loop indices
[2715]	367	INTEGER :: irgb, ierror, ioptio, nqsr ! local integer
[4147]	368	INTEGER :: ios ! Local integer output status for namelist read
[2715]	369	REAL(wp) :: zz0, zc0 , zc1, zcoef ! local scalars
	370	REAL(wp) :: zz1, zc2 , zc3, zchl ! - -
[3294]	371	REAL(wp), POINTER, DIMENSION(:,: ) :: zekb, zekg, zekr
	372	REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea
[2715]	373	!
[1423]	374	CHARACTER(len=100) :: cn_dir ! Root directory for location of ssr files
	375	TYPE(FLD_N) :: sn_chl ! informations about the chlorofyl field to be read
[2715]	376	!!
[4161]	377	NAMELIST/namtra_qsr/ sn_chl, cn_dir, ln_traqsr, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, ln_qsr_ice, &
[2528]	378	& nn_chldta, rn_abs, rn_si0, rn_si1
[3]	379	!!----------------------------------------------------------------------
	380
[3294]	381	!
	382	IF( nn_timing == 1 ) CALL timing_start('tra_qsr_init')
	383	!
[4161]	384	! clem init for oatte and iatte
[4205]	385	IF( .NOT. ln_rstart ) THEN
	386	oatte(:,:) = 1._wp
	387	iatte(:,:) = 1._wp
	388	ENDIF
[4161]	389	!
[3294]	390	CALL wrk_alloc( jpi, jpj, zekb, zekg, zekr )
	391	CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )
	392	!
[2715]	393
[4147]	394	REWIND( numnam_ref ) ! Namelist namtra_qsr in reference namelist : Ratio and length of penetration
	395	READ ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901)
	396	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_qsr in reference namelist', lwp )
	397
	398	REWIND( numnam_cfg ) ! Namelist namtra_qsr in configuration namelist : Ratio and length of penetration
	399	READ ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902 )
	400	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist', lwp )
	401	WRITE ( numond, namtra_qsr )
[1423]	402	!
	403	IF(lwp) THEN ! control print
	404	WRITE(numout,*)
	405	WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation'
	406	WRITE(numout,*) '~~~~~~~~~~~~'
[1601]	407	WRITE(numout,*) ' Namelist namtra_qsr : set the parameter of penetration'
	408	WRITE(numout,*) ' Light penetration (T) or not (F) ln_traqsr = ', ln_traqsr
	409	WRITE(numout,*) ' RGB (Red-Green-Blue) light penetration ln_qsr_rgb = ', ln_qsr_rgb
	410	WRITE(numout,*) ' 2 band light penetration ln_qsr_2bd = ', ln_qsr_2bd
	411	WRITE(numout,*) ' bio-model light penetration ln_qsr_bio = ', ln_qsr_bio
[4161]	412	WRITE(numout,*) ' light penetration for ice-model LIM3 ln_qsr_ice = ', ln_qsr_ice
[1601]	413	WRITE(numout,*) ' RGB : Chl data (=1) or cst value (=0) nn_chldta = ', nn_chldta
	414	WRITE(numout,*) ' RGB & 2 bands: fraction of light (rn_si1) rn_abs = ', rn_abs
	415	WRITE(numout,*) ' RGB & 2 bands: shortess depth of extinction rn_si0 = ', rn_si0
	416	WRITE(numout,*) ' 2 bands: longest depth of extinction rn_si1 = ', rn_si1
[4161]	417	WRITE(numout,*) ' light penetration for ice-model LIM3 ln_qsr_ice = ', ln_qsr_ice
[1423]	418	ENDIF
[1448]	419
	420	IF( ln_traqsr ) THEN ! control consistency
	421	!
[1601]	422	IF( .NOT.lk_qsr_bio .AND. ln_qsr_bio ) THEN
	423	CALL ctl_warn( 'No bio model : force ln_qsr_bio = FALSE ' )
[1448]	424	ln_qsr_bio = .FALSE.
	425	ENDIF
	426	!
	427	ioptio = 0 ! Parameter control
	428	IF( ln_qsr_rgb ) ioptio = ioptio + 1
	429	IF( ln_qsr_2bd ) ioptio = ioptio + 1
	430	IF( ln_qsr_bio ) ioptio = ioptio + 1
	431	!
[2528]	432	IF( ioptio /= 1 ) &
	433	CALL ctl_stop( ' Choose ONE type of light penetration in namelist namtra_qsr', &
	434	& ' 2 bands, 3 RGB bands or bio-model light penetration' )
[1448]	435	!
[1455]	436	IF( ln_qsr_rgb .AND. nn_chldta == 0 ) nqsr = 1
	437	IF( ln_qsr_rgb .AND. nn_chldta == 1 ) nqsr = 2
	438	IF( ln_qsr_2bd ) nqsr = 3
	439	IF( ln_qsr_bio ) nqsr = 4
	440	!
[1448]	441	IF(lwp) THEN ! Print the choice
	442	WRITE(numout,*)
[2528]	443	IF( nqsr == 1 ) WRITE(numout,*) ' R-G-B light penetration - Constant Chlorophyll'
	444	IF( nqsr == 2 ) WRITE(numout,*) ' R-G-B light penetration - Chl data '
	445	IF( nqsr == 3 ) WRITE(numout,*) ' 2 bands light penetration'
[1455]	446	IF( nqsr == 4 ) WRITE(numout,*) ' bio-model light penetration'
[1448]	447	ENDIF
	448	!
	449	ENDIF
[1423]	450	! ! ===================================== !
	451	IF( ln_traqsr ) THEN ! Initialisation of Light Penetration !
	452	! ! ===================================== !
	453	!
[2528]	454	xsi0r = 1.e0 / rn_si0
	455	xsi1r = 1.e0 / rn_si1
[1423]	456	! ! ---------------------------------- !
	457	IF( ln_qsr_rgb ) THEN ! Red-Green-Blue light penetration !
	458	! ! ---------------------------------- !
	459	!
[2528]	460	CALL trc_oce_rgb( rkrgb ) !* tabulated attenuation coef.
	461	!
	462	! !* level of light extinction
	463	IF( ln_sco ) THEN ; nksr = jpkm1
	464	ELSE ; nksr = trc_oce_ext_lev( r_si2, 0.33e2 )
[457]	465	ENDIF
[2528]	466
[4292]	467	IF(lwp) WRITE(numout,*) ' level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
[1423]	468	!
	469	IF( nn_chldta == 1 ) THEN !* Chl data : set sf_chl structure
	470	IF(lwp) WRITE(numout,*)
	471	IF(lwp) WRITE(numout,*) ' Chlorophyll read in a file'
	472	ALLOCATE( sf_chl(1), STAT=ierror )
	473	IF( ierror > 0 ) THEN
	474	CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_chl structure' ) ; RETURN
	475	ENDIF
[2528]	476	ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1) )
	477	IF( sn_chl%ln_tint )ALLOCATE( sf_chl(1)%fdta(jpi,jpj,1,2) )
[1423]	478	! ! fill sf_chl with sn_chl and control print
	479	CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init', &
[1601]	480	& 'Solar penetration function of read chlorophyll', 'namtra_qsr' )
[1423]	481	!
	482	ELSE !* constant Chl : compute once for all the distribution of light (etot3)
	483	IF(lwp) WRITE(numout,*)
	484	IF(lwp) WRITE(numout,*) ' Constant Chlorophyll concentration = 0.05'
[2528]	485	IF( lk_vvl ) THEN ! variable volume
	486	IF(lwp) WRITE(numout,*) ' key_vvl: light distribution will be computed at each time step'
	487	ELSE ! constant volume: computes one for all
	488	IF(lwp) WRITE(numout,*) ' fixed volume: light distribution computed one for all'
	489	!
	490	zchl = 0.05 ! constant chlorophyll
	491	irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
	492	zekb(:,:) = rkrgb(1,irgb) ! Separation in R-G-B depending of the chlorophyll
	493	zekg(:,:) = rkrgb(2,irgb)
	494	zekr(:,:) = rkrgb(3,irgb)
	495	!
	496	zcoef = ( 1. - rn_abs ) / 3.e0 ! equi-partition in R-G-B
	497	ze0(:,:,1) = rn_abs
	498	ze1(:,:,1) = zcoef
	499	ze2(:,:,1) = zcoef
	500	ze3(:,:,1) = zcoef
	501	zea(:,:,1) = tmask(:,:,1) ! = ( ze0+ze1+z2+ze3 ) * tmask
[1423]	502
[2528]	503	DO jk = 2, nksr+1
[1423]	504	!CDIR NOVERRCHK
[2528]	505	DO jj = 1, jpj
[1423]	506	!CDIR NOVERRCHK
[2528]	507	DO ji = 1, jpi
[4292]	508	zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * xsi0r )
	509	zc1 = ze1(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekb(ji,jj) )
	510	zc2 = ze2(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekg(ji,jj) )
	511	zc3 = ze3(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekr(ji,jj) )
[2528]	512	ze0(ji,jj,jk) = zc0
	513	ze1(ji,jj,jk) = zc1
	514	ze2(ji,jj,jk) = zc2
	515	ze3(ji,jj,jk) = zc3
	516	zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,jk)
	517	END DO
[1423]	518	END DO
[2528]	519	END DO
	520	!
	521	DO jk = 1, nksr
[3625]	522	etot3(:,:,jk) = r1_rau0_rcp * ( zea(:,:,jk) - zea(:,:,jk+1) )
[1423]	523	END DO
[2528]	524	etot3(:,:,nksr+1:jpk) = 0.e0 ! below 400m set to zero
	525	ENDIF
[1423]	526	ENDIF
	527	!
[1448]	528	ENDIF
[1423]	529	! ! ---------------------------------- !
[1448]	530	IF( ln_qsr_2bd ) THEN ! 2 bands light penetration !
[1423]	531	! ! ---------------------------------- !
	532	!
	533	! ! level of light extinction
	534	nksr = trc_oce_ext_lev( rn_si1, 1.e2 )
	535	IF(lwp) THEN
	536	WRITE(numout,*)
[4292]	537	IF(lwp) WRITE(numout,*) ' level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
[1423]	538	ENDIF
	539	!
[2528]	540	IF( lk_vvl ) THEN ! variable volume
	541	IF(lwp) WRITE(numout,*) ' key_vvl: light distribution will be computed at each time step'
	542	ELSE ! constant volume: computes one for all
[3625]	543	zz0 = rn_abs * r1_rau0_rcp
	544	zz1 = ( 1. - rn_abs ) * r1_rau0_rcp
[2528]	545	DO jk = 1, nksr !* solar heat absorbed at T-point computed once for all
	546	DO jj = 1, jpj ! top 400 meters
	547	DO ji = 1, jpi
	548	zc0 = zz0 * EXP( -fsdepw(ji,jj,jk )xsi0r ) + zz1 EXP( -fsdepw(ji,jj,jk )*xsi1r )
	549	zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)xsi0r ) + zz1 EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
	550	etot3(ji,jj,jk) = ( zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1) )
	551	END DO
[1423]	552	END DO
	553	END DO
[2528]	554	etot3(:,:,nksr+1:jpk) = 0.e0 ! below 400m set to zero
	555	!
	556	ENDIF
[187]	557	ENDIF
[1423]	558	! ! ===================================== !
	559	ELSE ! No light penetration !
	560	! ! ===================================== !
[457]	561	IF(lwp) THEN
	562	WRITE(numout,*)
	563	WRITE(numout,*) 'tra_qsr_init : NO solar flux penetration'
	564	WRITE(numout,*) '~~~~~~~~~~~~'
	565	ENDIF
[3]	566	ENDIF
[503]	567	!
[3294]	568	CALL wrk_dealloc( jpi, jpj, zekb, zekg, zekr )
	569	CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )
[2715]	570	!
[3294]	571	IF( nn_timing == 1 ) CALL timing_stop('tra_qsr_init')
	572	!
[3]	573	END SUBROUTINE tra_qsr_init
	574
	575	!!======================================================================
	576	END MODULE traqsr

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

Download in other formats: