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 NEMO/trunk/src/OCE/TRA – NEMO

Context Navigation

source: NEMO/trunk/src/OCE/TRA/traqsr.F90 @ 14045

Last change on this file since 14045 was 13982, checked in by smasson, 3 years ago
trunk: merge dev_r13923_Tiling_Cleanup_MPI3_LoopFusion into the trunk
Property svn:keywords set to `Id`
File size: 24.5 KB

Rev	Line
[3]	1	MODULE traqsr
	2	!!======================================================================
	3	!! * MODULE traqsr *
[6140]	4	!! Ocean physics: solar radiation penetration in the top ocean levels
[3]	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)
[6140]	12	!! 3.6 ! 2012-05 (C. Rousset) store attenuation coef for use in ice model
[6403]	13	!! 3.6 ! 2015-12 (O. Aumont, J. Jouanno, C. Ethe) use vertical profile of chlorophyll
[6140]	14	!! 3.7 ! 2015-11 (G. Madec, A. Coward) remove optimisation for fix volume
[3]	15	!!----------------------------------------------------------------------
[503]	16
	17	!!----------------------------------------------------------------------
[6140]	18	!! tra_qsr : temperature trend due to the penetration of solar radiation
	19	!! tra_qsr_init : initialization of the qsr penetration
[3]	20	!!----------------------------------------------------------------------
[6140]	21	USE oce ! ocean dynamics and active tracers
	22	USE phycst ! physical constants
	23	USE dom_oce ! ocean space and time domain
[13982]	24	USE domain, ONLY : dom_tile
[6140]	25	USE sbc_oce ! surface boundary condition: ocean
	26	USE trc_oce ! share SMS/Ocean variables
[4990]	27	USE trd_oce ! trends: ocean variables
	28	USE trdtra ! trends manager: tracers
[6140]	29	!
	30	USE in_out_manager ! I/O manager
	31	USE prtctl ! Print control
[9019]	32	USE iom ! I/O library
[6140]	33	USE fldread ! read input fields
	34	USE restart ! ocean restart
	35	USE lib_mpp ! MPP library
	36	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[3294]	37	USE timing ! Timing
[3]	38
	39	IMPLICIT NONE
	40	PRIVATE
	41
[2528]	42	PUBLIC tra_qsr ! routine called by step.F90 (ln_traqsr=T)
[5407]	43	PUBLIC tra_qsr_init ! routine called by nemogcm.F90
[3]	44
[4147]	45	! !!* Namelist namtra_qsr: penetrative solar radiation
	46	LOGICAL , PUBLIC :: ln_traqsr !: light absorption (qsr) flag
	47	LOGICAL , PUBLIC :: ln_qsr_rgb !: Red-Green-Blue light absorption flag
	48	LOGICAL , PUBLIC :: ln_qsr_2bd !: 2 band light absorption flag
	49	LOGICAL , PUBLIC :: ln_qsr_bio !: bio-model light absorption flag
	50	INTEGER , PUBLIC :: nn_chldta !: use Chlorophyll data (=1) or not (=0)
	51	REAL(wp), PUBLIC :: rn_abs !: fraction absorbed in the very near surface (RGB & 2 bands)
	52	REAL(wp), PUBLIC :: rn_si0 !: very near surface depth of extinction (RGB & 2 bands)
	53	REAL(wp), PUBLIC :: rn_si1 !: deepest depth of extinction (water type I) (2 bands)
[6140]	54	!
	55	INTEGER , PUBLIC :: nksr !: levels below which the light cannot penetrate (depth larger than 391 m)
[5407]	56
[6140]	57	INTEGER, PARAMETER :: np_RGB = 1 ! R-G-B light penetration with constant Chlorophyll
	58	INTEGER, PARAMETER :: np_RGBc = 2 ! R-G-B light penetration with Chlorophyll data
	59	INTEGER, PARAMETER :: np_2BD = 3 ! 2 bands light penetration
	60	INTEGER, PARAMETER :: np_BIO = 4 ! bio-model light penetration
	61	!
	62	INTEGER :: nqsr ! user choice of the type of light penetration
	63	REAL(wp) :: xsi0r ! inverse of rn_si0
	64	REAL(wp) :: xsi1r ! inverse of rn_si1
	65	!
[13333]	66	REAL(wp) , PUBLIC, DIMENSION(3,61) :: rkrgb ! tabulated attenuation coefficients for RGB absorption
[1423]	67	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_chl ! structure of input Chl (file informations, fields read)
[3]	68
	69	!! * Substitutions
[12377]	70	# include "do_loop_substitute.h90"
[13237]	71	# include "domzgr_substitute.h90"
[3]	72	!!----------------------------------------------------------------------
[9598]	73	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[888]	74	!! $Id$
[10068]	75	!! Software governed by the CeCILL license (see ./LICENSE)
[3]	76	!!----------------------------------------------------------------------
	77	CONTAINS
	78
[12377]	79	SUBROUTINE tra_qsr( kt, Kmm, pts, Krhs )
[3]	80	!!----------------------------------------------------------------------
	81	!! * ROUTINE tra_qsr *
	82	!!
	83	!! ** Purpose : Compute the temperature trend due to the solar radiation
[6140]	84	!! penetration and add it to the general temperature trend.
[3]	85	!!
[1423]	86	!! ** Method : The profile of the solar radiation within the ocean is defined
	87	!! through 2 wavebands (rn_si0,rn_si1) or 3 wavebands (RGB) and a ratio rn_abs
	88	!! Considering the 2 wavebands case:
	89	!! I(k) = Qsr( rn_absEXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) )
	90	!! The temperature trend associated with the solar radiation penetration
[12489]	91	!! is given by : zta = 1/e3t dk[ I ] / (rho0*Cp)
[3]	92	!! At the bottom, boudary condition for the radiation is no flux :
	93	!! all heat which has not been absorbed in the above levels is put
	94	!! in the last ocean level.
[6140]	95	!! The computation is only done down to the level where
	96	!! I(k) < 1.e-15 W/m2 (i.e. over the top nksr levels) .
[3]	97	!!
	98	!! ** Action : - update ta with the penetrative solar radiation trend
[6140]	99	!! - send trend for further diagnostics (l_trdtra=T)
[1423]	100	!!
	101	!! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
	102	!! Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516.
[6403]	103	!! Morel, A. et Berthon, JF, 1989, Limnol Oceanogr 34(8), 1545-1562
[503]	104	!!----------------------------------------------------------------------
[12377]	105	INTEGER, INTENT(in ) :: kt ! ocean time-step
	106	INTEGER, INTENT(in ) :: Kmm, Krhs ! time level indices
	107	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation
[2715]	108	!
[6140]	109	INTEGER :: ji, jj, jk ! dummy loop indices
[13982]	110	INTEGER :: irgb, isi, iei, isj, iej ! local integers
[6140]	111	REAL(wp) :: zchl, zcoef, z1_2 ! local scalars
	112	REAL(wp) :: zc0 , zc1 , zc2 , zc3 ! - -
[4161]	113	REAL(wp) :: zzc0, zzc1, zzc2, zzc3 ! - -
[13205]	114	REAL(wp) :: zz0 , zz1 , ze3t, zlui ! - -
	115	REAL(wp) :: zCb, zCmax, zpsi, zpsimax, zrdpsi, zCze
	116	REAL(wp) :: zlogc, zlogze, zlogCtot, zlogCze
	117	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ze0, ze1, ze2, ze3
	118	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt, zetot, ztmp3d
[3]	119	!!----------------------------------------------------------------------
[3294]	120	!
[9019]	121	IF( ln_timing ) CALL timing_start('tra_qsr')
[3294]	122	!
[13982]	123	IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile
	124	IF( kt == nit000 ) THEN
	125	IF(lwp) WRITE(numout,*)
	126	IF(lwp) WRITE(numout,*) 'tra_qsr : penetration of the surface solar radiation'
	127	IF(lwp) WRITE(numout,*) '~~~~~~~'
	128	ENDIF
[3]	129	ENDIF
[6140]	130	!
	131	IF( l_trdtra ) THEN ! trends diagnostic: save the input temperature trend
[13982]	132	ALLOCATE( ztrdt(jpi,jpj,jpk) )
[12377]	133	ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
[216]	134	ENDIF
[6140]	135	!
	136	! !-----------------------------------!
	137	! ! before qsr induced heat content !
	138	! !-----------------------------------!
[13982]	139	! NOTE: [tiling-comms-merge] Many DO loop bounds changed (probably more than necessary) to avoid changing results when using tiling. Some bounds were also adjusted to account for those changed in tra_atf
	140	IF( ntsi == Nis0 ) THEN ; isi = nn_hls ; ELSE ; isi = 0 ; ENDIF ! Avoid double-counting when using tiling
	141	IF( ntsj == Njs0 ) THEN ; isj = nn_hls ; ELSE ; isj = 0 ; ENDIF
	142	IF( ntei == Nie0 ) THEN ; iei = nn_hls ; ELSE ; iei = 0 ; ENDIF
	143	IF( ntej == Nje0 ) THEN ; iej = nn_hls ; ELSE ; iej = 0 ; ENDIF
	144
[6140]	145	IF( kt == nit000 ) THEN !== 1st time step ==!
[12489]	146	IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 .AND. .NOT.l_1st_euler ) THEN ! read in restart
[6140]	147	z1_2 = 0.5_wp
[13982]	148	IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile
	149	IF(lwp) WRITE(numout,*) ' nit000-1 qsr tracer content forcing field read in the restart file'
	150	CALL iom_get( numror, jpdom_auto, 'qsr_hc_b', qsr_hc_b ) ! before heat content trend due to Qsr flux
	151	ENDIF
[2528]	152	ELSE ! No restart or restart not found: Euler forward time stepping
[6140]	153	z1_2 = 1._wp
[13982]	154	DO_3D( isj, iej, isi, iei, 1, jpk )
	155	qsr_hc_b(ji,jj,jk) = 0._wp
	156	END_3D
[2528]	157	ENDIF
[6140]	158	ELSE !== Swap of qsr heat content ==!
	159	z1_2 = 0.5_wp
[13982]	160	DO_3D( isj, iej, isi, iei, 1, jpk )
	161	qsr_hc_b(ji,jj,jk) = qsr_hc(ji,jj,jk)
	162	END_3D
[2528]	163	ENDIF
[6140]	164	!
	165	! !--------------------------------!
	166	SELECT CASE( nqsr ) ! now qsr induced heat content !
	167	! !--------------------------------!
	168	!
	169	CASE( np_BIO ) !== bio-model fluxes ==!
	170	!
[13982]	171	DO_3D( isj, iej, isi, iei, 1, nksr )
	172	qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) )
	173	END_3D
[6140]	174	!
	175	CASE( np_RGB , np_RGBc ) !== R-G-B fluxes ==!
	176	!
[13982]	177	ALLOCATE( ze0 (A2D(nn_hls)) , ze1 (A2D(nn_hls)) , &
	178	& ze2 (A2D(nn_hls)) , ze3 (A2D(nn_hls)) , &
	179	& ztmp3d(A2D(nn_hls),nksr + 1) )
[6140]	180	!
	181	IF( nqsr == np_RGBc ) THEN !* Variable Chlorophyll
[13982]	182	IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only for the full domain
	183	IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 ) ! Use full domain
	184	CALL fld_read( kt, 1, sf_chl ) ! Read Chl data and provides it at the current time step
	185	IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 1 ) ! Revert to tile domain
	186	ENDIF
[13205]	187	!
	188	! Separation in R-G-B depending on the surface Chl
	189	! perform and store as many of the 2D calculations as possible
	190	! before the 3D loop (use the temporary 2D arrays to replace the
	191	! most expensive calculations)
	192	!
[13982]	193	DO_2D( isj, iej, isi, iei )
[13205]	194	! zlogc = log(zchl)
	195	zlogc = LOG ( MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) )
	196	! zc1 : log(zCze) = log (1.12 * zchl**0.803)
	197	zc1 = 0.113328685307 + 0.803 * zlogc
	198	! zc2 : log(zCtot) = log(40.6 * zchl**0.459)
	199	zc2 = 3.703768066608 + 0.459 * zlogc
	200	! zc3 : log(zze) = log(568.2 * zCtot**(-0.746))
	201	zc3 = 6.34247346942 - 0.746 * zc2
	202	! IF( log(zze) > log(102.) ) log(zze) = log(200.0 * zCtot**(-0.293))
	203	IF( zc3 > 4.62497281328 ) zc3 = 5.298317366548 - 0.293 * zc2
	204	!
	205	ze0(ji,jj) = zlogc ! ze0 = log(zchl)
	206	ze1(ji,jj) = EXP( zc1 ) ! ze1 = zCze
	207	ze2(ji,jj) = 1._wp / ( 0.710 + zlogc * ( 0.159 + zlogc * 0.021 ) ) ! ze2 = 1/zdelpsi
	208	ze3(ji,jj) = EXP( - zc3 ) ! ze3 = 1/zze
	209	END_2D
	210
	211	!
[13982]	212	DO_3D( isj, iej, isi, iei, 1, nksr + 1 )
[13205]	213	! zchl = ALOG( ze0(ji,jj) )
	214	zlogc = ze0(ji,jj)
	215	!
	216	zCb = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) )
	217	zCmax = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 )
	218	zpsimax = 0.6 - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) )
	219	! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 )
	220	!
	221	zCze = ze1(ji,jj)
	222	zrdpsi = ze2(ji,jj) ! 1/zdelpsi
	223	zpsi = ze3(ji,jj) * gdepw(ji,jj,jk,Kmm) ! gdepw/zze
	224	!
	225	! NB. make sure zchl value is such that: zchl = MIN( 10. , MAX( 0.03, zchl ) )
	226	zchl = MIN( 10. , MAX( 0.03, zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) ) )
	227	! Convert chlorophyll value to attenuation coefficient look-up table index
	228	ztmp3d(ji,jj,jk) = 41 + 20.*LOG10(zchl) + 1.e-15
	229	END_3D
	230	ELSE !* constant chlorophyll
	231	zchl = 0.05
	232	! NB. make sure constant value is such that:
	233	zchl = MIN( 10. , MAX( 0.03, zchl ) )
	234	! Convert chlorophyll value to attenuation coefficient look-up table index
	235	zlui = 41 + 20.*LOG10(zchl) + 1.e-15
[6403]	236	DO jk = 1, nksr + 1
[13982]	237	ztmp3d(:,:,jk) = zlui
[3]	238	END DO
[4161]	239	ENDIF
[1423]	240	!
[6140]	241	zcoef = ( 1. - rn_abs ) / 3._wp !* surface equi-partition in R-G-B
[13982]	242	DO_2D( isj, iej, isi, iei )
[13205]	243	ze0(ji,jj) = rn_abs * qsr(ji,jj)
	244	ze1(ji,jj) = zcoef * qsr(ji,jj)
	245	ze2(ji,jj) = zcoef * qsr(ji,jj)
	246	ze3(ji,jj) = zcoef * qsr(ji,jj)
	247	! store the surface SW radiation; re-use the surface ztmp3d array
	248	! since the surface attenuation coefficient is not used
	249	ztmp3d(ji,jj,1) = qsr(ji,jj)
[12377]	250	END_2D
[6140]	251	!
[13497]	252	! !* interior equi-partition in R-G-B depending on vertical profile of Chl
[13982]	253	DO_3D( isj, iej, isi, iei, 2, nksr + 1 )
[13205]	254	ze3t = e3t(ji,jj,jk-1,Kmm)
	255	irgb = NINT( ztmp3d(ji,jj,jk) )
	256	zc0 = ze0(ji,jj) * EXP( - ze3t * xsi0r )
	257	zc1 = ze1(ji,jj) * EXP( - ze3t * rkrgb(1,irgb) )
	258	zc2 = ze2(ji,jj) * EXP( - ze3t * rkrgb(2,irgb) )
	259	zc3 = ze3(ji,jj) * EXP( - ze3t * rkrgb(3,irgb) )
	260	ze0(ji,jj) = zc0
	261	ze1(ji,jj) = zc1
	262	ze2(ji,jj) = zc2
	263	ze3(ji,jj) = zc3
	264	ztmp3d(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk)
	265	END_3D
[6140]	266	!
[13982]	267	DO_3D( isj, iej, isi, iei, 1, nksr ) !* now qsr induced heat content
[13205]	268	qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( ztmp3d(ji,jj,jk) - ztmp3d(ji,jj,jk+1) )
[12377]	269	END_3D
[187]	270	!
[13205]	271	DEALLOCATE( ze0 , ze1 , ze2 , ze3 , ztmp3d )
[6140]	272	!
	273	CASE( np_2BD ) !== 2-bands fluxes ==!
	274	!
[12489]	275	zz0 = rn_abs * r1_rho0_rcp ! surface equi-partition in 2-bands
	276	zz1 = ( 1. - rn_abs ) * r1_rho0_rcp
[13982]	277	DO_3D( isj, iej, isi, iei, 1, nksr ) !* now qsr induced heat content
[12377]	278	zc0 = zz0 * EXP( -gdepw(ji,jj,jk ,Kmm)xsi0r ) + zz1 EXP( -gdepw(ji,jj,jk ,Kmm)*xsi1r )
	279	zc1 = zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)xsi0r ) + zz1 EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi1r )
	280	qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) )
	281	END_3D
[2528]	282	!
[6140]	283	END SELECT
	284	!
	285	! !-----------------------------!
[13497]	286	! ! update to the temp. trend !
	287	! !-----------------------------!
[13295]	288	DO_3D( 0, 0, 0, 0, 1, nksr )
[12377]	289	pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) &
[13237]	290	& + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) &
	291	& / e3t(ji,jj,jk,Kmm)
[12377]	292	END_3D
[6140]	293	!
[9019]	294	! sea-ice: store the 1st ocean level attenuation coefficient
[13982]	295	DO_2D( isj, iej, isi, iei )
[12489]	296	IF( qsr(ji,jj) /= 0._wp ) THEN ; fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) )
[12377]	297	ELSE ; fraqsr_1lev(ji,jj) = 1._wp
	298	ENDIF
	299	END_2D
[2528]	300	!
[13982]	301	! TEMP: [tiling] This change not necessary and working array can use A2D(nn_hls) if using XIOS (subdomain support)
	302	IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only for the full domain
	303	IF( iom_use('qsr3d') ) THEN ! output the shortwave Radiation distribution
	304	ALLOCATE( zetot(jpi,jpj,jpk) )
	305	zetot(:,:,nksr+1:jpk) = 0._wp ! below ~400m set to zero
	306	DO jk = nksr, 1, -1
	307	zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * rho0_rcp
	308	END DO
	309	CALL iom_put( 'qsr3d', zetot ) ! 3D distribution of shortwave Radiation
	310	DEALLOCATE( zetot )
	311	ENDIF
[2528]	312	ENDIF
[6140]	313	!
[13982]	314	IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only on the last tile
	315	IF( lrst_oce ) THEN ! write in the ocean restart file
	316	CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b' , qsr_hc )
	317	CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev )
	318	ENDIF
[6140]	319	ENDIF
	320	!
[503]	321	IF( l_trdtra ) THEN ! qsr tracers trends saved for diagnostics
[12377]	322	ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
	323	CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt )
[13982]	324	DEALLOCATE( ztrdt )
[3]	325	ENDIF
[457]	326	! ! print mean trends (used for debugging)
[12377]	327	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' qsr - Ta: ', mask1=tmask, clinfo3='tra-ta' )
[503]	328	!
[9019]	329	IF( ln_timing ) CALL timing_stop('tra_qsr')
[3294]	330	!
[3]	331	END SUBROUTINE tra_qsr
	332
	333
	334	SUBROUTINE tra_qsr_init
	335	!!----------------------------------------------------------------------
	336	!! * ROUTINE tra_qsr_init *
	337	!!
	338	!! ** Purpose : Initialization for the penetrative solar radiation
	339	!!
	340	!! ** Method : The profile of solar radiation within the ocean is set
[1423]	341	!! from two length scale of penetration (rn_si0,rn_si1) and a ratio
[1601]	342	!! (rn_abs). These parameters are read in the namtra_qsr namelist. The
[3]	343	!! default values correspond to clear water (type I in Jerlov'
	344	!! (1968) classification.
	345	!! called by tra_qsr at the first timestep (nit000)
	346	!!
[1423]	347	!! ** Action : - initialize rn_si0, rn_si1 and rn_abs
[3]	348	!!
[503]	349	!! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
[3]	350	!!----------------------------------------------------------------------
[6140]	351	INTEGER :: ji, jj, jk ! dummy loop indices
	352	INTEGER :: ios, irgb, ierror, ioptio ! local integer
	353	REAL(wp) :: zz0, zc0 , zc1, zcoef ! local scalars
	354	REAL(wp) :: zz1, zc2 , zc3, zchl ! - -
[2715]	355	!
[1423]	356	CHARACTER(len=100) :: cn_dir ! Root directory for location of ssr files
	357	TYPE(FLD_N) :: sn_chl ! informations about the chlorofyl field to be read
[2715]	358	!!
[9019]	359	NAMELIST/namtra_qsr/ sn_chl, cn_dir, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, &
[2528]	360	& nn_chldta, rn_abs, rn_si0, rn_si1
[3]	361	!!----------------------------------------------------------------------
[3294]	362	!
[6140]	363	READ ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901)
[11536]	364	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_qsr in reference namelist' )
[3294]	365	!
[4147]	366	READ ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902 )
[11536]	367	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist' )
[4624]	368	IF(lwm) WRITE ( numond, namtra_qsr )
[1423]	369	!
	370	IF(lwp) THEN ! control print
	371	WRITE(numout,*)
	372	WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation'
	373	WRITE(numout,*) '~~~~~~~~~~~~'
[1601]	374	WRITE(numout,*) ' Namelist namtra_qsr : set the parameter of penetration'
[6140]	375	WRITE(numout,*) ' RGB (Red-Green-Blue) light penetration ln_qsr_rgb = ', ln_qsr_rgb
	376	WRITE(numout,*) ' 2 band light penetration ln_qsr_2bd = ', ln_qsr_2bd
	377	WRITE(numout,*) ' bio-model light penetration ln_qsr_bio = ', ln_qsr_bio
	378	WRITE(numout,*) ' RGB : Chl data (=1) or cst value (=0) nn_chldta = ', nn_chldta
	379	WRITE(numout,*) ' RGB & 2 bands: fraction of light (rn_si1) rn_abs = ', rn_abs
	380	WRITE(numout,*) ' RGB & 2 bands: shortess depth of extinction rn_si0 = ', rn_si0
	381	WRITE(numout,*) ' 2 bands: longest depth of extinction rn_si1 = ', rn_si1
	382	WRITE(numout,*)
[1423]	383	ENDIF
[6140]	384	!
	385	ioptio = 0 ! Parameter control
	386	IF( ln_qsr_rgb ) ioptio = ioptio + 1
	387	IF( ln_qsr_2bd ) ioptio = ioptio + 1
	388	IF( ln_qsr_bio ) ioptio = ioptio + 1
	389	!
	390	IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE type of light penetration in namelist namtra_qsr', &
	391	& ' 2 bands, 3 RGB bands or bio-model light penetration' )
	392	!
	393	IF( ln_qsr_rgb .AND. nn_chldta == 0 ) nqsr = np_RGB
	394	IF( ln_qsr_rgb .AND. nn_chldta == 1 ) nqsr = np_RGBc
	395	IF( ln_qsr_2bd ) nqsr = np_2BD
	396	IF( ln_qsr_bio ) nqsr = np_BIO
	397	!
	398	! ! Initialisation
	399	xsi0r = 1._wp / rn_si0
	400	xsi1r = 1._wp / rn_si1
	401	!
	402	SELECT CASE( nqsr )
	403	!
	404	CASE( np_RGB , np_RGBc ) !== Red-Green-Blue light penetration ==!
	405	!
[9169]	406	IF(lwp) WRITE(numout,*) ' ==>>> R-G-B light penetration '
[6140]	407	!
	408	CALL trc_oce_rgb( rkrgb ) ! tabulated attenuation coef.
	409	!
	410	nksr = trc_oce_ext_lev( r_si2, 33._wp ) ! level of light extinction
	411	!
	412	IF(lwp) WRITE(numout,*) ' level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
	413	!
	414	IF( nqsr == np_RGBc ) THEN ! Chl data : set sf_chl structure
[9169]	415	IF(lwp) WRITE(numout,*) ' ==>>> Chlorophyll read in a file'
[6140]	416	ALLOCATE( sf_chl(1), STAT=ierror )
	417	IF( ierror > 0 ) THEN
	418	CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_chl structure' ) ; RETURN
	419	ENDIF
	420	ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1) )
	421	IF( sn_chl%ln_tint ) ALLOCATE( sf_chl(1)%fdta(jpi,jpj,1,2) )
	422	! ! fill sf_chl with sn_chl and control print
	423	CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init', &
[7646]	424	& 'Solar penetration function of read chlorophyll', 'namtra_qsr' , no_print )
[1448]	425	ENDIF
[6140]	426	IF( nqsr == np_RGB ) THEN ! constant Chl
[9169]	427	IF(lwp) WRITE(numout,*) ' ==>>> Constant Chlorophyll concentration = 0.05'
[6140]	428	ENDIF
[1448]	429	!
[6140]	430	CASE( np_2BD ) !== 2 bands light penetration ==!
[1448]	431	!
[9169]	432	IF(lwp) WRITE(numout,*) ' ==>>> 2 bands light penetration'
[1448]	433	!
[6140]	434	nksr = trc_oce_ext_lev( rn_si1, 100._wp ) ! level of light extinction
	435	IF(lwp) WRITE(numout,*) ' level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
[1455]	436	!
[6140]	437	CASE( np_BIO ) !== BIO light penetration ==!
[1448]	438	!
[9169]	439	IF(lwp) WRITE(numout,*) ' ==>>> bio-model light penetration'
[7646]	440	IF( .NOT.lk_top ) CALL ctl_stop( 'No bio model : ln_qsr_bio = true impossible ' )
[1423]	441	!
[13333]	442	CALL trc_oce_rgb( rkrgb ) ! tabulated attenuation coef.
	443	!
	444	nksr = trc_oce_ext_lev( r_si2, 33._wp ) ! level of light extinction
	445	!
	446	IF(lwp) WRITE(numout,*) ' level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
	447	!
[6140]	448	END SELECT
[503]	449	!
[7753]	450	qsr_hc(:,:,:) = 0._wp ! now qsr heat content set to zero where it will not be computed
[6140]	451	!
	452	! 1st ocean level attenuation coefficient (used in sbcssm)
[5407]	453	IF( iom_varid( numror, 'fraqsr_1lev', ldstop = .FALSE. ) > 0 ) THEN
[13970]	454	CALL iom_get( numror, jpdom_auto, 'fraqsr_1lev' , fraqsr_1lev )
[5407]	455	ELSE
[7753]	456	fraqsr_1lev(:,:) = 1._wp ! default : no penetration
[5407]	457	ENDIF
	458	!
[3]	459	END SUBROUTINE tra_qsr_init
	460
	461	!!======================================================================
	462	END MODULE traqsr

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

Download in other formats: