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_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90 @ 4332

Last change on this file since 4332 was 4332, checked in by clem, 10 years ago
update LIM3 to fix remaining bugs. Now working in global and regional config.
Property svn:keywords set to `Id`
File size: 30.0 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)
[4038]	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
[4038]	30	USE restart ! ocean restart
[2715]	31	USE lib_mpp ! MPP library
[3294]	32	USE wrk_nemo ! Memory Allocation
	33	USE timing ! Timing
[4220]	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
[1601]	42	! !!* Namelist namtra_qsr: penetrative solar radiation
[1423]	43	LOGICAL , PUBLIC :: ln_traqsr = .TRUE. !: light absorption (qsr) flag
	44	LOGICAL , PUBLIC :: ln_qsr_rgb = .FALSE. !: Red-Green-Blue light absorption flag
[1448]	45	LOGICAL , PUBLIC :: ln_qsr_2bd = .TRUE. !: 2 band light absorption flag
[1423]	46	LOGICAL , PUBLIC :: ln_qsr_bio = .FALSE. !: bio-model light absorption flag
[4038]	47	LOGICAL , PUBLIC :: ln_qsr_ice = .FALSE. !: light penetration for ice-model LIM3 (clem)
[1423]	48	INTEGER , PUBLIC :: nn_chldta = 0 !: use Chlorophyll data (=1) or not (=0)
	49	REAL(wp), PUBLIC :: rn_abs = 0.58_wp !: fraction absorbed in the very near surface (RGB & 2 bands)
	50	REAL(wp), PUBLIC :: rn_si0 = 0.35_wp !: very near surface depth of extinction (RGB & 2 bands)
	51	REAL(wp), PUBLIC :: rn_si1 = 23.0_wp !: 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 ! - -
[4038]	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
[4038]	164	! clem: store attenuation coefficient of the first ocean level
[4220]	165	IF ( lk_lim3 .AND. ln_qsr_ice ) THEN
[4038]	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
[4038]	233	! clem: store attenuation coefficient of the first ocean level
[4220]	234	IF ( lk_lim3 .AND. ln_qsr_ice ) THEN
[4038]	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
[4038]	257	! clem: store attenuation coefficient of the first ocean level
[4220]	258	IF ( lk_lim3 .AND. ln_qsr_ice ) THEN
[4038]	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
[4038]	281	! clem: store attenuation coefficient of the first ocean level
[4220]	282	IF ( lk_lim3 .AND. ln_qsr_ice ) THEN
[4038]	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
[4038]	300	! clem: store attenuation coefficient of the first ocean level
[4220]	301	IF ( lk_lim3 .AND. ln_qsr_ice ) THEN
[4038]	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	!
[4220]	366	INTEGER :: ji, jj, jk ! dummy loop indices
[2715]	367	INTEGER :: irgb, ierror, ioptio, nqsr ! local integer
	368	REAL(wp) :: zz0, zc0 , zc1, zcoef ! local scalars
	369	REAL(wp) :: zz1, zc2 , zc3, zchl ! - -
[3294]	370	REAL(wp), POINTER, DIMENSION(:,: ) :: zekb, zekg, zekr
	371	REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea
[2715]	372	!
[1423]	373	CHARACTER(len=100) :: cn_dir ! Root directory for location of ssr files
	374	TYPE(FLD_N) :: sn_chl ! informations about the chlorofyl field to be read
[2715]	375	!!
[4038]	376	NAMELIST/namtra_qsr/ sn_chl, cn_dir, ln_traqsr, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, ln_qsr_ice, &
[2528]	377	& nn_chldta, rn_abs, rn_si0, rn_si1
[3]	378	!!----------------------------------------------------------------------
	379
[3294]	380	!
	381	IF( nn_timing == 1 ) CALL timing_start('tra_qsr_init')
	382	!
[4038]	383	! clem init for oatte and iatte
[4220]	384	IF( .NOT. ln_rstart ) THEN
	385	oatte(:,:) = 1._wp
	386	iatte(:,:) = 1._wp
	387	ENDIF
[4038]	388	!
[3294]	389	CALL wrk_alloc( jpi, jpj, zekb, zekg, zekr )
	390	CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )
	391	!
[2715]	392
[1423]	393	cn_dir = './' ! directory in which the model is executed
	394	! ... default values (NB: frequency positive => hours, negative => months)
	395	! ! file ! frequency ! variable ! time interp ! clim ! 'yearly' or ! weights ! rotation !
	396	! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! filename ! pairs !
	397	sn_chl = FLD_N( 'chlorophyll' , -1 , 'CHLA' , .true. , .true. , 'yearly' , '' , '' )
	398	!
[1601]	399	REWIND( numnam ) ! Read Namelist namtra_qsr : ratio and length of penetration
	400	READ ( numnam, namtra_qsr )
[1423]	401	!
	402	IF(lwp) THEN ! control print
	403	WRITE(numout,*)
	404	WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation'
	405	WRITE(numout,*) '~~~~~~~~~~~~'
[1601]	406	WRITE(numout,*) ' Namelist namtra_qsr : set the parameter of penetration'
	407	WRITE(numout,*) ' Light penetration (T) or not (F) ln_traqsr = ', ln_traqsr
	408	WRITE(numout,*) ' RGB (Red-Green-Blue) light penetration ln_qsr_rgb = ', ln_qsr_rgb
	409	WRITE(numout,*) ' 2 band light penetration ln_qsr_2bd = ', ln_qsr_2bd
	410	WRITE(numout,*) ' bio-model light penetration ln_qsr_bio = ', ln_qsr_bio
[4038]	411	WRITE(numout,*) ' light penetration for ice-model LIM3 ln_qsr_ice = ', ln_qsr_ice
[1601]	412	WRITE(numout,*) ' RGB : Chl data (=1) or cst value (=0) nn_chldta = ', nn_chldta
	413	WRITE(numout,*) ' RGB & 2 bands: fraction of light (rn_si1) rn_abs = ', rn_abs
	414	WRITE(numout,*) ' RGB & 2 bands: shortess depth of extinction rn_si0 = ', rn_si0
	415	WRITE(numout,*) ' 2 bands: longest depth of extinction rn_si1 = ', rn_si1
[4220]	416	WRITE(numout,*) ' light penetration for ice-model LIM3 ln_qsr_ice = ', ln_qsr_ice
[1423]	417	ENDIF
[1448]	418
	419	IF( ln_traqsr ) THEN ! control consistency
	420	!
[1601]	421	IF( .NOT.lk_qsr_bio .AND. ln_qsr_bio ) THEN
	422	CALL ctl_warn( 'No bio model : force ln_qsr_bio = FALSE ' )
[1448]	423	ln_qsr_bio = .FALSE.
	424	ENDIF
	425	!
	426	ioptio = 0 ! Parameter control
	427	IF( ln_qsr_rgb ) ioptio = ioptio + 1
	428	IF( ln_qsr_2bd ) ioptio = ioptio + 1
	429	IF( ln_qsr_bio ) ioptio = ioptio + 1
	430	!
[2528]	431	IF( ioptio /= 1 ) &
	432	CALL ctl_stop( ' Choose ONE type of light penetration in namelist namtra_qsr', &
	433	& ' 2 bands, 3 RGB bands or bio-model light penetration' )
[1448]	434	!
[1455]	435	IF( ln_qsr_rgb .AND. nn_chldta == 0 ) nqsr = 1
	436	IF( ln_qsr_rgb .AND. nn_chldta == 1 ) nqsr = 2
	437	IF( ln_qsr_2bd ) nqsr = 3
	438	IF( ln_qsr_bio ) nqsr = 4
	439	!
[1448]	440	IF(lwp) THEN ! Print the choice
	441	WRITE(numout,*)
[2528]	442	IF( nqsr == 1 ) WRITE(numout,*) ' R-G-B light penetration - Constant Chlorophyll'
	443	IF( nqsr == 2 ) WRITE(numout,*) ' R-G-B light penetration - Chl data '
	444	IF( nqsr == 3 ) WRITE(numout,*) ' 2 bands light penetration'
[1455]	445	IF( nqsr == 4 ) WRITE(numout,*) ' bio-model light penetration'
[1448]	446	ENDIF
	447	!
	448	ENDIF
[1423]	449	! ! ===================================== !
	450	IF( ln_traqsr ) THEN ! Initialisation of Light Penetration !
	451	! ! ===================================== !
	452	!
[2528]	453	xsi0r = 1.e0 / rn_si0
	454	xsi1r = 1.e0 / rn_si1
[1423]	455	! ! ---------------------------------- !
	456	IF( ln_qsr_rgb ) THEN ! Red-Green-Blue light penetration !
	457	! ! ---------------------------------- !
	458	!
[2528]	459	CALL trc_oce_rgb( rkrgb ) !* tabulated attenuation coef.
	460	!
	461	! !* level of light extinction
	462	IF( ln_sco ) THEN ; nksr = jpkm1
	463	ELSE ; nksr = trc_oce_ext_lev( r_si2, 0.33e2 )
[457]	464	ENDIF
[2528]	465
	466	IF(lwp) WRITE(numout,*) ' level of light extinction = ', nksr, ' ref depth = ', gdepw_0(nksr+1), ' m'
[1423]	467	!
	468	IF( nn_chldta == 1 ) THEN !* Chl data : set sf_chl structure
	469	IF(lwp) WRITE(numout,*)
	470	IF(lwp) WRITE(numout,*) ' Chlorophyll read in a file'
	471	ALLOCATE( sf_chl(1), STAT=ierror )
	472	IF( ierror > 0 ) THEN
	473	CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_chl structure' ) ; RETURN
	474	ENDIF
[2528]	475	ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1) )
	476	IF( sn_chl%ln_tint )ALLOCATE( sf_chl(1)%fdta(jpi,jpj,1,2) )
[1423]	477	! ! fill sf_chl with sn_chl and control print
	478	CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init', &
[1601]	479	& 'Solar penetration function of read chlorophyll', 'namtra_qsr' )
[1423]	480	!
	481	ELSE !* constant Chl : compute once for all the distribution of light (etot3)
	482	IF(lwp) WRITE(numout,*)
	483	IF(lwp) WRITE(numout,*) ' Constant Chlorophyll concentration = 0.05'
[2528]	484	IF( lk_vvl ) THEN ! variable volume
	485	IF(lwp) WRITE(numout,*) ' key_vvl: light distribution will be computed at each time step'
	486	ELSE ! constant volume: computes one for all
	487	IF(lwp) WRITE(numout,*) ' fixed volume: light distribution computed one for all'
	488	!
	489	zchl = 0.05 ! constant chlorophyll
	490	irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
	491	zekb(:,:) = rkrgb(1,irgb) ! Separation in R-G-B depending of the chlorophyll
	492	zekg(:,:) = rkrgb(2,irgb)
	493	zekr(:,:) = rkrgb(3,irgb)
	494	!
	495	zcoef = ( 1. - rn_abs ) / 3.e0 ! equi-partition in R-G-B
	496	ze0(:,:,1) = rn_abs
	497	ze1(:,:,1) = zcoef
	498	ze2(:,:,1) = zcoef
	499	ze3(:,:,1) = zcoef
	500	zea(:,:,1) = tmask(:,:,1) ! = ( ze0+ze1+z2+ze3 ) * tmask
[1423]	501
[2528]	502	DO jk = 2, nksr+1
[1423]	503	!CDIR NOVERRCHK
[2528]	504	DO jj = 1, jpj
[1423]	505	!CDIR NOVERRCHK
[2528]	506	DO ji = 1, jpi
	507	zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t_0(ji,jj,jk-1) * xsi0r )
	508	zc1 = ze1(ji,jj,jk-1) * EXP( - fse3t_0(ji,jj,jk-1) * zekb(ji,jj) )
	509	zc2 = ze2(ji,jj,jk-1) * EXP( - fse3t_0(ji,jj,jk-1) * zekg(ji,jj) )
	510	zc3 = ze3(ji,jj,jk-1) * EXP( - fse3t_0(ji,jj,jk-1) * zekr(ji,jj) )
	511	ze0(ji,jj,jk) = zc0
	512	ze1(ji,jj,jk) = zc1
	513	ze2(ji,jj,jk) = zc2
	514	ze3(ji,jj,jk) = zc3
	515	zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,jk)
	516	END DO
[1423]	517	END DO
[2528]	518	END DO
	519	!
	520	DO jk = 1, nksr
[3625]	521	etot3(:,:,jk) = r1_rau0_rcp * ( zea(:,:,jk) - zea(:,:,jk+1) )
[1423]	522	END DO
[2528]	523	etot3(:,:,nksr+1:jpk) = 0.e0 ! below 400m set to zero
	524	ENDIF
[1423]	525	ENDIF
	526	!
[1448]	527	ENDIF
[1423]	528	! ! ---------------------------------- !
[1448]	529	IF( ln_qsr_2bd ) THEN ! 2 bands light penetration !
[1423]	530	! ! ---------------------------------- !
	531	!
	532	! ! level of light extinction
	533	nksr = trc_oce_ext_lev( rn_si1, 1.e2 )
	534	IF(lwp) THEN
	535	WRITE(numout,*)
[2528]	536	IF(lwp) WRITE(numout,*) ' level of light extinction = ', nksr, ' ref depth = ', gdepw_0(nksr+1), ' m'
[1423]	537	ENDIF
	538	!
[2528]	539	IF( lk_vvl ) THEN ! variable volume
	540	IF(lwp) WRITE(numout,*) ' key_vvl: light distribution will be computed at each time step'
	541	ELSE ! constant volume: computes one for all
[3625]	542	zz0 = rn_abs * r1_rau0_rcp
	543	zz1 = ( 1. - rn_abs ) * r1_rau0_rcp
[2528]	544	DO jk = 1, nksr !* solar heat absorbed at T-point computed once for all
	545	DO jj = 1, jpj ! top 400 meters
	546	DO ji = 1, jpi
	547	zc0 = zz0 * EXP( -fsdepw(ji,jj,jk )xsi0r ) + zz1 EXP( -fsdepw(ji,jj,jk )*xsi1r )
	548	zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)xsi0r ) + zz1 EXP( -fsdepw(ji,jj,jk+1)*xsi1r )
	549	etot3(ji,jj,jk) = ( zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1) )
	550	END DO
[1423]	551	END DO
	552	END DO
[2528]	553	etot3(:,:,nksr+1:jpk) = 0.e0 ! below 400m set to zero
	554	!
	555	ENDIF
[187]	556	ENDIF
[1423]	557	! ! ===================================== !
	558	ELSE ! No light penetration !
	559	! ! ===================================== !
[457]	560	IF(lwp) THEN
	561	WRITE(numout,*)
	562	WRITE(numout,*) 'tra_qsr_init : NO solar flux penetration'
	563	WRITE(numout,*) '~~~~~~~~~~~~'
	564	ENDIF
[3]	565	ENDIF
[503]	566	!
[3294]	567	CALL wrk_dealloc( jpi, jpj, zekb, zekg, zekr )
	568	CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )
[2715]	569	!
[3294]	570	IF( nn_timing == 1 ) CALL timing_stop('tra_qsr_init')
	571	!
[3]	572	END SUBROUTINE tra_qsr_init
	573
	574	!!======================================================================
	575	END MODULE traqsr

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

Download in other formats: