Changeset 14990
- Timestamp:
- 2021-06-14T21:19:43+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/TRA/traqsr.F90
r14618 r14990 13 13 !! 3.6 ! 2015-12 (O. Aumont, J. Jouanno, C. Ethe) use vertical profile of chlorophyll 14 14 !! 3.7 ! 2015-11 (G. Madec, A. Coward) remove optimisation for fix volume 15 !! 4.0 ! 2020-11 (A. Coward) optimisation 16 !! 4.5 ! 2021-03 (G. Madec) further optimisation + adaptation for RK3 15 17 !!---------------------------------------------------------------------- 16 18 17 19 !!---------------------------------------------------------------------- 18 20 !! tra_qsr : temperature trend due to the penetration of solar radiation 21 !! qsr_RGBc : IR + RGB light penetration with Chlorophyll data case 22 !! qsr_RGB : IR + RGB light penetration with constant Chlorophyll case 23 !! qsr_2BD : 2 bands (InfraRed + Visible light) case 24 !! qsr_ext_lev : level of extinction for each bands 19 25 !! tra_qsr_init : initialization of the qsr penetration 20 26 !!---------------------------------------------------------------------- … … 53 59 REAL(wp), PUBLIC :: rn_si1 !: deepest depth of extinction (water type I) (2 bands) 54 60 ! 55 INTEGER , PUBLIC :: nksr !: levels below which the light cannot penetrate (depth larger than 391 m)56 57 61 INTEGER, PARAMETER :: np_RGB = 1 ! R-G-B light penetration with constant Chlorophyll 58 62 INTEGER, PARAMETER :: np_RGBc = 2 ! R-G-B light penetration with Chlorophyll data … … 60 64 INTEGER, PARAMETER :: np_BIO = 4 ! bio-model light penetration 61 65 ! 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 66 INTEGER :: nqsr ! user choice of the type of light penetration 67 INTEGER :: nc_rgb ! RGB with cst Chlorophyll: index associated with the chosen Chl value 68 ! 69 ! ! extinction level 70 INTEGER :: nk0 !: IR (depth larger ~12 m) 71 INTEGER :: nkV !: Visible light (depth larger than ~840 m) 72 INTEGER :: nkR, nkG, nkB !: RGB (depth larger than ~100 m, ~470 m, ~1700 m, resp.) 73 ! 74 INTEGER, PUBLIC :: nksr !: =nkV, i.e. maximum level of light extinction (used in traatf(_qco).F90) 75 ! 76 ! ! inverse of attenuation length 77 REAL(wp) :: r1_si0 ! all schemes : infrared = 1/rn_si0 78 REAL(wp) :: r1_si1 ! 2 band : mean RGB = 1/rn_si1 79 REAL(wp) :: r1_LR, r1_LG, r1_LB ! RGB with constant Chl (np_RGB) 65 80 ! 66 81 REAL(wp) , PUBLIC, DIMENSION(3,61) :: rkrgb ! tabulated attenuation coefficients for RGB absorption … … 77 92 CONTAINS 78 93 79 SUBROUTINE tra_qsr( kt, Kmm, pts, Krhs , kstg)94 SUBROUTINE tra_qsr( kt, Kmm, pts, Krhs ) 80 95 !!---------------------------------------------------------------------- 81 96 !! *** ROUTINE tra_qsr *** … … 85 100 !! 86 101 !! ** 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_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) ) 90 !! The temperature trend associated with the solar radiation penetration 91 !! is given by : zta = 1/e3t dk[ I ] / (rho0*Cp) 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. 102 !! through 2 wavebands (rn_si0,rn_si1) or 3 wavebands (RGB) or computed by 103 !! the biogeochemical model 95 104 !! The computation is only done down to the level where 96 !! I(k) < 1.e-15 W/m2 (i.e. over the top nk srlevels) .97 !! 98 !! ** Action : - update t awith the penetrative solar radiation trend105 !! I(k) < 1.e-15 W/m2 (i.e. over the top nk levels) . 106 !! 107 !! ** Action : - update ts(jp_tem) with the penetrative solar radiation trend 99 108 !! - send trend for further diagnostics (l_trdtra=T) 100 !! 101 !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp. 102 !! Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516. 103 !! Morel, A. et Berthon, JF, 1989, Limnol Oceanogr 34(8), 1545-1562 104 !!---------------------------------------------------------------------- 105 INTEGER, INTENT(in ) :: kt, Kmm, Krhs ! ocean time-step and time level indices 109 !!---------------------------------------------------------------------- 110 INTEGER, INTENT(in ) :: kt, Kmm, Krhs ! ocean time-step and time-level indices 106 111 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation 107 INTEGER , OPTIONAL , INTENT(in ) :: kstg ! RK3 stage index108 112 ! 109 113 INTEGER :: ji, jj, jk ! dummy loop indices 110 INTEGER :: irgb, isi, iei, isj, iej ! local integers 111 INTEGER :: istg_1, istg_3 ! - - 112 REAL(wp) :: zchl, zcoef, z1_2 ! local scalars 113 REAL(wp) :: zc0 , zc1 , zc2 , zc3 ! - - 114 REAL(wp) :: zzc0, zzc1, zzc2, zzc3 ! - - 115 REAL(wp) :: zz0 , zz1 , ze3t, zlui ! - - 116 REAL(wp) :: zCb, zCmax, zpsi, zpsimax, zrdpsi, zCze 117 REAL(wp) :: zlogc, zlogze, zlogCtot, zlogCze 118 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ze0, ze1, ze2, ze3 119 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt, zetot, ztmp3d 114 REAL(wp) :: z1_2, ze3t ! local scalars 115 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt, zetot 120 116 !!---------------------------------------------------------------------- 121 117 ! 122 118 IF( ln_timing ) CALL timing_start('tra_qsr') 123 !124 IF( PRESENT( kstg ) ) THEN ! RK3 : a few things have to be done at only a specific stage125 istg_1 = kstg ; istg_3 = kstg126 ELSE ! MLF : only one call by time step127 istg_1 = 1 ; istg_3 = 3128 ENDIF129 119 ! 130 120 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile … … 136 126 ENDIF 137 127 ! 138 IF( l_trdtra ) THEN 128 IF( l_trdtra ) THEN ! trends diagnostic: save the input temperature trend 139 129 ALLOCATE( ztrdt(jpi,jpj,jpk) ) 140 130 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) 141 131 ENDIF 142 ! 143 ! !-----------------------------------! 144 ! ! before qsr induced heat content ! 145 ! !-----------------------------------! 146 IF( ntsi == Nis0 ) THEN ; isi = nn_hls ; ELSE ; isi = 0 ; ENDIF ! Avoid double-counting when using tiling 147 IF( ntsj == Njs0 ) THEN ; isj = nn_hls ; ELSE ; isj = 0 ; ENDIF 148 IF( ntei == Nie0 ) THEN ; iei = nn_hls ; ELSE ; iei = 0 ; ENDIF 149 IF( ntej == Nje0 ) THEN ; iej = nn_hls ; ELSE ; iej = 0 ; ENDIF 150 151 IF( kt == nit000 .AND. istg_1 == 1 ) THEN !== 1st time step ==! (RK3: only at stage 1) 132 ! 133 #if ! defined key_RK3 134 ! ! MLF only : heat content trend due to Qsr flux (qsr_hc) 135 ! 136 IF( kt == nit000 ) THEN !== 1st time step ==! 152 137 IF( ln_rstart .AND. .NOT.l_1st_euler ) THEN ! read in restart 153 138 z1_2 = 0.5_wp … … 158 143 ELSE ! No restart or Euler forward at 1st time step 159 144 z1_2 = 1._wp 160 DO_3D( isi, iei, isj, iej, 1, jpk )145 DO_3D( 0,0, 0,0, 1, jpk ) 161 146 qsr_hc_b(ji,jj,jk) = 0._wp 162 147 END_3D 163 148 ENDIF 164 ELSE IF( istg_3 == 3 ) THEN!== Swap of qsr heat content ==!149 ELSE !== Swap of qsr heat content ==! 165 150 z1_2 = 0.5_wp 166 DO_3D( isi, iei, isj, iej, 1, jpk )151 DO_3D( 0,0, 0,0, 1, jpk ) 167 152 qsr_hc_b(ji,jj,jk) = qsr_hc(ji,jj,jk) 168 153 END_3D 169 154 ENDIF 170 ! 171 ! !--------------------------------! 172 SELECT CASE( nqsr ) ! now qsr induced heat content ! 173 ! !--------------------------------! 174 ! 175 CASE( np_BIO ) !== bio-model fluxes ==! 176 ! 177 DO_3D( isi, iei, isj, iej, 1, nksr ) 155 #endif 156 157 ! !----------------------------! 158 SELECT CASE( nqsr ) ! qsr induced heat content ! 159 ! !----------------------------! 160 ! 161 CASE( np_RGBc ) ; CALL qsr_RGBc( kt, Kmm, pts, Krhs ) !== R-G-B fluxes using chlorophyll data ==! with Morel &Berthon (1989) vertical profile 162 ! 163 CASE( np_RGB ) ; CALL qsr_RGB ( kt, Kmm, pts, Krhs ) !== R-G-B fluxes with constant chlorophyll ==! 164 ! 165 CASE( np_2BD ) ; CALL qsr_2BD ( Kmm, pts, Krhs ) !== 2-bands fluxes ==! 166 ! 167 CASE( np_BIO ) !== bio-model fluxes ==! 168 DO_3D( 0, 0, 0, 0, 1, nkV ) 169 #if defined key_RK3 170 ! !- RK3 : temperature trend at jk t-level 171 ze3t = e3t(ji,jj,jk,Kmm) 172 pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) / ze3t 173 #else 174 ! !- MLF : heat content trend due to Qsr flux (qsr_hc) 178 175 qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) 176 #endif 179 177 END_3D 180 ! 181 CASE( np_RGB , np_RGBc ) !== R-G-B fluxes ==! 182 ! 183 ALLOCATE( ze0 (A2D(nn_hls)) , ze1 (A2D(nn_hls)) , & 184 & ze2 (A2D(nn_hls)) , ze3 (A2D(nn_hls)) , & 185 & ztmp3d(A2D(nn_hls),nksr + 1) ) 186 ! 187 IF( nqsr == np_RGBc ) THEN !* Variable Chlorophyll 188 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only for the full domain 189 IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 ) ! Use full domain 190 CALL fld_read( kt, 1, sf_chl ) ! Read Chl data and provides it at the current time step 191 IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 1 ) ! Revert to tile domain 192 ENDIF 193 ! 194 ! Separation in R-G-B depending on the surface Chl 195 ! perform and store as many of the 2D calculations as possible 196 ! before the 3D loop (use the temporary 2D arrays to replace the 197 ! most expensive calculations) 198 ! 199 DO_2D( isi, iei, isj, iej ) 200 ! zlogc = log(zchl) 201 zlogc = LOG ( MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) ) 202 ! zc1 : log(zCze) = log (1.12 * zchl**0.803) 203 zc1 = 0.113328685307 + 0.803 * zlogc 204 ! zc2 : log(zCtot) = log(40.6 * zchl**0.459) 205 zc2 = 3.703768066608 + 0.459 * zlogc 206 ! zc3 : log(zze) = log(568.2 * zCtot**(-0.746)) 207 zc3 = 6.34247346942 - 0.746 * zc2 208 ! IF( log(zze) > log(102.) ) log(zze) = log(200.0 * zCtot**(-0.293)) 209 IF( zc3 > 4.62497281328 ) zc3 = 5.298317366548 - 0.293 * zc2 210 ! 211 ze0(ji,jj) = zlogc ! ze0 = log(zchl) 212 ze1(ji,jj) = EXP( zc1 ) ! ze1 = zCze 213 ze2(ji,jj) = 1._wp / ( 0.710 + zlogc * ( 0.159 + zlogc * 0.021 ) ) ! ze2 = 1/zdelpsi 214 ze3(ji,jj) = EXP( - zc3 ) ! ze3 = 1/zze 215 END_2D 216 ! 217 DO_3D( isi, iei, isj, iej, 1, nksr + 1 ) 218 ! zchl = ALOG( ze0(ji,jj) ) 219 zlogc = ze0(ji,jj) 220 ! 221 zCb = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) ) 222 zCmax = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 ) 223 zpsimax = 0.6 - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) ) 224 ! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 ) 225 ! 226 zCze = ze1(ji,jj) 227 zrdpsi = ze2(ji,jj) ! 1/zdelpsi 228 zpsi = ze3(ji,jj) * gdepw(ji,jj,jk,Kmm) ! gdepw/zze 229 ! 230 ! NB. make sure zchl value is such that: zchl = MIN( 10. , MAX( 0.03, zchl ) ) 231 zchl = MIN( 10. , MAX( 0.03, zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) ) ) 232 ! Convert chlorophyll value to attenuation coefficient look-up table index 233 ztmp3d(ji,jj,jk) = 41 + 20.*LOG10(zchl) + 1.e-15 234 END_3D 235 ELSE !* constant chlorophyll 236 zchl = 0.05 237 ! NB. make sure constant value is such that: 238 zchl = MIN( 10. , MAX( 0.03, zchl ) ) 239 ! Convert chlorophyll value to attenuation coefficient look-up table index 240 zlui = 41 + 20.*LOG10(zchl) + 1.e-15 241 DO jk = 1, nksr + 1 242 ztmp3d(:,:,jk) = zlui 178 ! !- sea-ice : store the 1st level attenuation coefficient 179 WHERE( etot3(A2D(0),1) /= 0._wp ) ; fraqsr_1lev(A2D(0)) = 1._wp - etot3(A2D(0),2) / etot3(A2D(0),1) 180 ELSEWHERE ; fraqsr_1lev(A2D(0)) = 1._wp 181 END WHERE 182 ! 183 END SELECT 184 ! 185 #if defined key_RK3 186 ! ! RK3 : diagnostics/output 187 IF( l_trdtra .OR. iom_use('qsr3d') ) THEN ! qsr diagnostics 188 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) 189 ! ! qsr tracers trends saved for diagnostics 190 IF( l_trdtra ) CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt ) 191 IF( iom_use('qsr3d') ) THEN ! qsr distribution 192 DO jk = nkV, 1, -1 193 ztrdt(:,:,jk) = ztrdt(:,:,jk+1) + qsr_hc(:,:,jk) * rho0_rcp 243 194 END DO 195 CALL iom_put( 'qsr3d', ztrdt ) ! 3D distribution of shortwave Radiation 244 196 ENDIF 245 ! 246 zcoef = ( 1. - rn_abs ) / 3._wp !* surface equi-partition in R-G-B 247 DO_2D( isi, iei, isj, iej ) 248 ze0(ji,jj) = rn_abs * qsr(ji,jj) 249 ze1(ji,jj) = zcoef * qsr(ji,jj) 250 ze2(ji,jj) = zcoef * qsr(ji,jj) 251 ze3(ji,jj) = zcoef * qsr(ji,jj) 252 ! store the surface SW radiation; re-use the surface ztmp3d array 253 ! since the surface attenuation coefficient is not used 254 ztmp3d(ji,jj,1) = qsr(ji,jj) 255 END_2D 256 ! 257 ! !* interior equi-partition in R-G-B depending on vertical profile of Chl 258 DO_3D( isi, iei, isj, iej, 2, nksr + 1 ) 259 ze3t = e3t(ji,jj,jk-1,Kmm) 260 irgb = NINT( ztmp3d(ji,jj,jk) ) 261 zc0 = ze0(ji,jj) * EXP( - ze3t * xsi0r ) 262 zc1 = ze1(ji,jj) * EXP( - ze3t * rkrgb(1,irgb) ) 263 zc2 = ze2(ji,jj) * EXP( - ze3t * rkrgb(2,irgb) ) 264 zc3 = ze3(ji,jj) * EXP( - ze3t * rkrgb(3,irgb) ) 265 ze0(ji,jj) = zc0 266 ze1(ji,jj) = zc1 267 ze2(ji,jj) = zc2 268 ze3(ji,jj) = zc3 269 ztmp3d(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk) 270 END_3D 271 ! 272 DO_3D( isi, iei, isj, iej, 1, nksr ) !* now qsr induced heat content 273 qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( ztmp3d(ji,jj,jk) - ztmp3d(ji,jj,jk+1) ) 274 END_3D 275 ! 276 DEALLOCATE( ze0 , ze1 , ze2 , ze3 , ztmp3d ) 277 ! 278 CASE( np_2BD ) !== 2-bands fluxes ==! 279 ! 280 zz0 = rn_abs * r1_rho0_rcp ! surface equi-partition in 2-bands 281 zz1 = ( 1. - rn_abs ) * r1_rho0_rcp 282 DO_3D( isi, iei, isj, iej, 1, nksr ) !* now qsr induced heat content 283 zc0 = zz0 * EXP( -gdepw(ji,jj,jk ,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk ,Kmm)*xsi1r ) 284 zc1 = zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi1r ) 285 qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) ) 286 END_3D 287 ! 288 END SELECT 289 ! 290 ! !-----------------------------! 291 ! ! update to the temp. trend ! 292 ! !-----------------------------! 197 DEALLOCATE( ztrdt ) 198 ENDIF 199 #else 200 ! ! MLF : add the temperature trend 293 201 DO_3D( 0, 0, 0, 0, 1, nksr ) 294 202 pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) & … … 297 205 END_3D 298 206 ! 207 !!st7-2 299 208 ! sea-ice: store the 1st ocean level attenuation coefficient 300 DO_2D( isi, iei, isj, iej)209 DO_2D( 0, 0, 0, 0 ) 301 210 IF( qsr(ji,jj) /= 0._wp ) THEN ; fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) ) 302 211 ELSE ; fraqsr_1lev(ji,jj) = 1._wp 303 212 ENDIF 304 213 END_2D 214 !!st7-2 end 305 215 ! 306 216 ! TEMP: [tiling] This change not necessary and working array can use A2D(nn_hls) if using XIOS (subdomain support) … … 317 227 ENDIF 318 228 ! 229 IF( l_trdtra ) THEN ! qsr tracers trends saved for diagnostics 230 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) 231 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt ) 232 DEALLOCATE( ztrdt ) 233 ENDIF 234 #endif 235 ! 319 236 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only on the last tile 320 237 IF( lrst_oce ) THEN ! write in the ocean restart file … … 324 241 ENDIF 325 242 ! 326 IF( l_trdtra ) THEN ! qsr tracers trends saved for diagnostics327 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)328 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt )329 DEALLOCATE( ztrdt )330 ENDIF243 !!st IF( l_trdtra ) THEN ! qsr tracers trends saved for diagnostics 244 !!st ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) 245 !!st CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt ) 246 !!st DEALLOCATE( ztrdt ) 247 !!st ENDIF 331 248 ! ! print mean trends (used for debugging) 332 249 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' qsr - Ta: ', mask1=tmask, clinfo3='tra-ta' ) … … 335 252 ! 336 253 END SUBROUTINE tra_qsr 254 255 256 SUBROUTINE qsr_RGBc( kt, Kmm, pts, Krhs ) 257 !!---------------------------------------------------------------------- 258 !! *** ROUTINE qsr_RGBc *** 259 !! 260 !! ** Purpose : Red-Green-Blue solar radiation using chlorophyll data 261 !! 262 !! ** Method : The profile of the solar radiation within the ocean is defined 263 !! through 2 wavebands (rn_si0,rn_si1) or 3 wavebands (RGB) and a ratio rn_abs 264 !! Considering the 2 wavebands case: 265 !! I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) ) 266 !! The temperature trend associated with the solar radiation penetration 267 !! is given by : zta = 1/e3t dk[ I ] / (rho0*Cp) 268 !! At the bottom, boudary condition for the radiation is no flux : 269 !! all heat which has not been absorbed in the above levels is put 270 !! in the last ocean level. 271 !! The computation is only done down to the level where 272 !! I(k) < 1.e-15 W/m2 (i.e. over the top nk levels) . 273 !! 274 !! ** Action : - update ta with the penetrative solar radiation trend 275 !! - send trend for further diagnostics (l_trdtra=T) 276 !! 277 !! Reference : Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516. 278 !! Morel, A. et Berthon, JF, 1989, Limnol Oceanogr 34(8), 1545-1562 279 !!---------------------------------------------------------------------- 280 INTEGER, INTENT(in ) :: kt, Kmm, Krhs ! ocean time-step and time-level indices 281 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation 282 !! 283 INTEGER :: ji, jj, jk ! dummy loop indices 284 INTEGER :: irgb ! local integer 285 REAL(wp) :: zc1 , zc2 , zc3, zchl ! local scalars 286 REAL(wp) :: zze0, zzeR, zzeG, zzeB, zzeT ! - - 287 REAL(wp) :: zz0 , zz1 , ze3t ! - - 288 REAL(wp) :: zCb, zCmax, zpsi, zpsimax, zrdpsi, zCze ! - - 289 REAL(wp) :: zlogc, zlogze, zlogCtot, zlogCze ! - - 290 REAL(wp), DIMENSION(A2D(0) ) :: ze0, zeR, zeG, zeB, zeT 291 REAL(wp), DIMENSION(A2D(0),0:3) :: zc 292 !!---------------------------------------------------------------------- 293 ! 294 ! 295 ! !===========================================! 296 ! !== R-G-B fluxes using chlorophyll data ==! with Morel &Berthon (1989) vertical profile 297 ! !===================================****====! 298 ! 299 ! != Chlorophyll data =! 300 ! 301 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only for the full domain 302 IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 ) ! Use full domain 303 CALL fld_read( kt, 1, sf_chl ) ! Read Chl data and provides it at the current time step 304 IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 1 ) ! Revert to tile domain 305 ENDIF 306 ! 307 DO_2D( 0, 0, 0, 0 ) ! pre-calculated expensive coefficient 308 zlogc = LOG( MAX( 0.03_wp, MIN( sf_chl(1)%fnow(ji,jj,1) ,10._wp ) ) ) ! zlogc = log(zchl) with 0.03 <= Chl >= 10. 309 zc1 = 0.113328685307 + 0.803 * zlogc ! zc1 : log(zCze) = log (1.12 * zchl**0.803) 310 zc2 = 3.703768066608 + 0.459 * zlogc ! zc2 : log(zCtot) = log(40.6 * zchl**0.459) 311 zc3 = 6.34247346942 - 0.746 * zc2 ! zc3 : log(zze) = log(568.2 * zCtot**(-0.746)) 312 IF( zc3 > 4.62497281328 ) zc3 = 5.298317366548 - 0.293 * zc2 ! IF(log(zze)>log(102)) log(zze) = log(200*zCtot**(-0.293)) 313 ! 314 zc(ji,jj,0) = zlogc ! ze(0) = log(zchl) 315 zc(ji,jj,1) = EXP( zc1 ) ! ze(1) = zCze 316 zc(ji,jj,2) = 1._wp / ( 0.710 + zlogc * ( 0.159 + zlogc * 0.021 ) ) ! ze(2) = 1/zdelpsi 317 zc(ji,jj,3) = EXP( - zc3 ) ! ze(3) = 1/zze 318 END_2D 319 ! 320 ! != surface light =! 321 ! 322 zz0 = rn_abs ! Infrared absorption 323 zz1 = ( 1._wp - rn_abs ) / 3._wp ! R-G-B equi-partition 324 ! 325 DO_2D( 0, 0, 0, 0 ) ! surface light 326 ze0(ji,jj) = zz0 * qsr(ji,jj) ; zeR(ji,jj) = zz1 * qsr(ji,jj) ! IR ; Red 327 zeG(ji,jj) = zz1 * qsr(ji,jj) ; zeB(ji,jj) = zz1 * qsr(ji,jj) ! Green ; Blue 328 zeT(ji,jj) = qsr(ji,jj) ! Total 329 END_2D 330 ! 331 ! != interior light =! 332 ! 333 DO jk = 1, nk0 !* near surface layers *! (< ~12 meters : IR + RGB ) 334 DO_2D( 0, 0, 0, 0 ) 335 ! !- inverse of RGB attenuation lengths 336 zlogc = zc(ji,jj,0) 337 zCb = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) ) 338 zCmax = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 ) 339 zpsimax = 0.6 - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) ) 340 ! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 ) 341 zCze = zc(ji,jj,1) 342 zrdpsi = zc(ji,jj,2) ! 1/zdelpsi 343 !!st05 zpsi = zc(ji,jj,3) * gdepw(ji,jj,jk,Kmm) ! gdepw/zze 344 zpsi = zc(ji,jj,3) * gdepw(ji,jj,jk+1,Kmm) ! gdepw/zze 345 ! ! make sure zchl value is such that: 0.03 < zchl < 10. 346 zchl = MAX( 0.03_wp , MIN( zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) , 10._wp ) ) 347 ! ! Convert chlorophyll value to attenuation coefficient 348 irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) ! look-up table index 349 ! Red ! Green ! Blue 350 r1_LR = rkrgb(3,irgb) ; r1_LG = rkrgb(2,irgb) ; r1_LB = rkrgb(1,irgb) 351 ! 352 ! !- fluxes at jk+1 w-level 353 ze3t = e3t(ji,jj,jk,Kmm) 354 zze0 = ze0(ji,jj) * EXP( - ze3t*r1_si0 ) ; zzeR = zeR(ji,jj) * EXP( - ze3t*r1_LR ) ! IR ; Red at jk+1 w-level 355 zzeG = zeG(ji,jj) * EXP( - ze3t*r1_LG ) ; zzeB = zeB(ji,jj) * EXP( - ze3t*r1_LB ) ! Green ; Blue - - 356 zzeT = ( zze0 + zzeB + zzeG + zzeR ) * wmask(ji,jj,jk+1) ! Total - - 357 !!st01 zzeT = ( zze0 + zzeR + zzeG + zzeB ) * wmask(ji,jj,jk+1) ! Total - - 358 ! 359 #if defined key_RK3 360 ! !- RK3 : temperature trend at jk t-level 361 pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t 362 #else 363 ! !- MLF : heat content trend due to Qsr flux (qsr_hc) 364 qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) 365 #endif 366 ze0(ji,jj) = zze0 ; zeR(ji,jj) = zzeR ! IR ; Red store at jk+1 w-level 367 zeG(ji,jj) = zzeG ; zeB(ji,jj) = zzeB ! Green ; Blue - - - 368 zeT(ji,jj) = zzeT ! total - - - 369 END_2D 370 !!stbug IF( jk == 1 ) THEN !- sea-ice : store the 1st level attenuation coefficient 371 !!stbug WHERE( qsr(A2D(0)) /= 0._wp ) ; fraqsr_1lev(A2D(0)) = 1._wp - zeT(A2D(0)) / qsr(A2D(0)) 372 !!stbug ELSEWHERE ; fraqsr_1lev(A2D(0)) = 1._wp 373 !!stbug END WHERE 374 !!stbug ENDIF 375 END DO 376 ! 377 IF( kt == nit000 ) THEN 378 IF(lwp) WRITE(numout,*) 'nk0+1= ', nk0+1, ' qsr IR max = ' , MAXVAL(ze0(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(ze0(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 379 IF(lwp) WRITE(numout,*) ' ', nk0+1, ' qsr R max = ' , MAXVAL(zeR(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeR(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 380 IF(lwp) WRITE(numout,*) ' ', nk0+1, ' qsr G max = ' , MAXVAL(zeG(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeG(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 381 IF(lwp) WRITE(numout,*) ' ', nk0+1, ' qsr B max = ' , MAXVAL(zeB(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeB(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 382 IF(lwp) WRITE(numout,*) ' ', nk0+1, ' qsr T max = ' , MAXVAL(zeT(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeT(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 383 ENDIF 384 ! 385 DO jk = nk0+1, nkR !* down to Red extinction *! (< ~71 meters : RGB , IR removed from calculation) 386 DO_2D( 0, 0, 0, 0 ) 387 ! !- inverse of RGB attenuation lengths 388 zlogc = zc(ji,jj,0) 389 zCb = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) ) 390 zCmax = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 ) 391 zpsimax = 0.6 - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) ) 392 ! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 ) 393 zCze = zc(ji,jj,1) 394 zrdpsi = zc(ji,jj,2) ! 1/zdelpsi 395 zpsi = zc(ji,jj,3) * gdepw(ji,jj,jk+1,Kmm) ! gdepw/zze 396 !!st05 zpsi = zc(ji,jj,3) * gdepw(ji,jj,jk,Kmm) ! gdepw/zze 397 ! ! make sure zchl value is such that: 0.03 < zchl < 10. 398 zchl = MAX( 0.03_wp , MIN( zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) , 10._wp ) ) 399 ! ! Convert chlorophyll value to attenuation coefficient 400 irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) ! look-up table index 401 ! Red ! Green ! Blue 402 r1_LR = rkrgb(3,irgb) ; r1_LG = rkrgb(2,irgb) ; r1_LB = rkrgb(1,irgb) 403 ! 404 ! !- fluxes at jk+1 w-level 405 ze3t = e3t(ji,jj,jk,Kmm) 406 zzeR = zeR(ji,jj) * EXP( - ze3t*r1_LR ) ! Red at jk+1 w-level 407 zzeG = zeG(ji,jj) * EXP( - ze3t*r1_LG ) ; zzeB = zeB(ji,jj) * EXP( - ze3t*r1_LB ) ! Green ; Blue - - 408 zzeT = ( zzeR + zzeG + zzeB ) * wmask(ji,jj,jk+1) ! Total - - 409 ! 410 #if defined key_RK3 411 ! !- RK3 : temperature trend at jk t-level 412 pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t 413 #else 414 ! !- MLF : heat content trend due to Qsr flux (qsr_hc) 415 qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) 416 #endif 417 zeR(ji,jj) = zzeR ! Red store at jk+1 w-level 418 zeG(ji,jj) = zzeG ; zeB(ji,jj) = zzeB ! Green ; Blue - - - 419 zeT(ji,jj) = zzeT ! total - - - 420 END_2D 421 END DO 422 ! 423 IF( kt == nit000 ) THEN 424 IF(lwp) WRITE(numout,*) 'nkR+1= ', nkR+1, ' qsr R max = ' , MAXVAL(zeR(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeR(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 425 IF(lwp) WRITE(numout,*) ' ', nkR+1, ' qsr G max = ' , MAXVAL(zeG(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeG(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 426 IF(lwp) WRITE(numout,*) ' ', nkR+1, ' qsr B max = ' , MAXVAL(zeB(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeB(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 427 IF(lwp) WRITE(numout,*) ' ', nkR+1, ' qsr T max = ' , MAXVAL(zeT(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeT(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 428 ENDIF 429 ! 430 DO jk = nkR+1, nkG !* down to Green extinction *! (< ~350 m : GB , IR+R removed from calculation) 431 DO_2D( 0, 0, 0, 0 ) 432 ! !- inverse of RGB attenuation lengths 433 zlogc = zc(ji,jj,0) 434 zCb = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) ) 435 zCmax = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 ) 436 zpsimax = 0.6 - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) ) 437 ! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 ) 438 zCze = zc(ji,jj,1) 439 zrdpsi = zc(ji,jj,2) ! 1/zdelpsi 440 zpsi = zc(ji,jj,3) * gdepw(ji,jj,jk+1,Kmm) ! gdepw/zze 441 !!st05 zpsi = zc(ji,jj,3) * gdepw(ji,jj,jk,Kmm) ! gdepw/zze 442 ! ! make sure zchl value is such that: 0.03 < zchl < 10. 443 zchl = MAX( 0.03_wp , MIN( zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) , 10._wp ) ) 444 ! ! Convert chlorophyll value to attenuation coefficient 445 irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) ! look-up table index 446 ! Green ! Blue 447 r1_LG = rkrgb(2,irgb) ; r1_LB = rkrgb(1,irgb) 448 ! 449 ! !- fluxes at jk+1 w-level 450 ze3t = e3t(ji,jj,jk,Kmm) 451 zzeG = zeG(ji,jj) * EXP( - ze3t * r1_LG ) ; zzeB = zeB(ji,jj) * EXP( - ze3t * r1_LB ) ! Green ; Blue 452 zzeT = ( zzeG + zzeB ) * wmask(ji,jj,jk+1) ! Total - - 453 #if defined key_RK3 454 ! !- RK3 : temperature trend at jk t-level 455 pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t 456 #else 457 ! !- MLF : heat content trend due to Qsr flux (qsr_hc) 458 qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) 459 #endif 460 zeG(ji,jj) = zzeG ; zeB(ji,jj) = zzeB ! Green ; Blue store at jk+1 w-level 461 zeT(ji,jj) = zzeT ! total - - - 462 END_2D 463 END DO 464 ! 465 IF( kt == nit000 ) THEN 466 IF(lwp) WRITE(numout,*) 'nkG+1= ', nkG+1, ' qsr G max = ' , MAXVAL(zeG(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeG(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 467 IF(lwp) WRITE(numout,*) ' ', nkG+1, ' qsr B max = ' , MAXVAL(zeB(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeB(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 468 IF(lwp) WRITE(numout,*) ' ', nkG+1, ' qsr T max = ' , MAXVAL(zeT(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeT(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 469 ENDIF 470 ! 471 DO jk = nkG+1, nkB !* down to Blue extinction *! (< ~1300 m : B , IR+RG removed from calculation) 472 DO_2D( 0, 0, 0, 0 ) 473 ! !- inverse of RGB attenuation lengths 474 zlogc = zc(ji,jj,0) 475 zCb = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) ) 476 zCmax = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 ) 477 zpsimax = 0.6 - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) ) 478 ! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 ) 479 zCze = zc(ji,jj,1) 480 zrdpsi = zc(ji,jj,2) ! 1/zdelpsi 481 zpsi = zc(ji,jj,3) * gdepw(ji,jj,jk+1,Kmm) ! gdepw/zze 482 !!st05 zpsi = zc(ji,jj,3) * gdepw(ji,jj,jk,Kmm) ! gdepw/zze 483 ! ! make sure zchl value is such that: 0.03 < zchl < 10. 484 zchl = MAX( 0.03_wp , MIN( zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) , 10._wp ) ) 485 ! ! Convert chlorophyll value to attenuation coefficient 486 irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) ! look-up table index 487 r1_LB = rkrgb(1,irgb) ! Blue 488 ! 489 ! !- fluxes at jk+1 w-level 490 ze3t = e3t(ji,jj,jk,Kmm) 491 zzeB = zeB(ji,jj) * EXP( - ze3t * r1_LB ) ! Blue 492 zzeT = ( zzeB ) * wmask(ji,jj,jk+1) ! Total - - 493 #if defined key_RK3 494 ! !- RK3 : temperature trend at jk t-level 495 pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t 496 #else 497 ! !- MLF : heat content trend due to Qsr flux (qsr_hc) 498 qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) 499 #endif 500 zeB(ji,jj) = zzeB ! Blue store at jk+1 w-level 501 zeT(ji,jj) = zzeT ! total - - - 502 END_2D 503 END DO 504 ! 505 IF( kt == nit000 ) THEN 506 IF(lwp) WRITE(numout,*) 'nkB+1= ', nkB+1, ' qsr T max = ' , MAXVAL(zeT), ' W/m2' , MAXVAL(zeT(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)), ' K/s' 507 ENDIF 508 ! 509 END SUBROUTINE qsr_RGBc 510 511 512 SUBROUTINE qsr_RGB( kt, Kmm, pts, Krhs ) 513 !!---------------------------------------------------------------------- 514 !! *** ROUTINE qsr_RGB *** 515 !! 516 !! ** Purpose : Red-Green-Blue solar radiation with constant chlorophyll 517 !! 518 !! ** Method : The profile of the solar radiation within the ocean is defined 519 !! through 2 wavebands (rn_si0,rn_si1) or 1 (rn_si0,rn_abs) + 3 wavebands (RGB) 520 !! At the bottom, boudary condition for the radiation is no flux : 521 !! all heat which has not been absorbed in the above levels is put 522 !! in the last ocean level. 523 !! For each band, the computation is only done down to the level where 524 !! I(k) < 1.e-15 W/m2 (i.e. over the top nk levels) . 525 !! 526 !! ** Action : - update ta with the penetrative solar radiation trend 527 !! - send trend for further diagnostics (l_trdtra=T) 528 !! 529 !! Reference : Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516. 530 !! Morel, A. et Berthon, JF, 1989, Limnol Oceanogr 34(8), 1545-1562 531 !!---------------------------------------------------------------------- 532 INTEGER, INTENT(in ) :: kt, Kmm, Krhs ! ocean time-step and time-level indices 533 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation 534 !! 535 INTEGER :: ji, jj, jk ! dummy loop indices 536 REAL(wp) :: zze0, zzeR, zzeG, zzeB, zzeT ! - - 537 REAL(wp) :: zz0 , zz1 , ze3t ! - - 538 REAL(wp), DIMENSION(A2D(0)) :: ze0, zeR, zeG, zeB, zeT 539 !!---------------------------------------------------------------------- 540 ! 541 ! 542 ! !==============================================! 543 ! !== R-G-B fluxes with constant chlorophyll ==! 544 ! !======================********================! 545 ! 546 ! != surface light =! 547 ! 548 zz0 = rn_abs ! Infrared absorption 549 zz1 = ( 1._wp - rn_abs ) / 3._wp ! surface equi-partition in R-G-B 550 ! 551 DO_2D( 0, 0, 0, 0 ) ! surface light 552 ze0(ji,jj) = zz0 * qsr(ji,jj) ; zeR(ji,jj) = zz1 * qsr(ji,jj) ! IR ; Red 553 zeG(ji,jj) = zz1 * qsr(ji,jj) ; zeB(ji,jj) = zz1 * qsr(ji,jj) ! Green ; Blue 554 zeT(ji,jj) = qsr(ji,jj) ! Total 555 END_2D 556 ! 557 ! != interior light =! 558 ! 559 DO jk = 1, nk0 !* near surface layers *! (< ~12 meters : IR + RGB ) 560 DO_2D( 0, 0, 0, 0 ) 561 ze3t = e3t(ji,jj,jk,Kmm) 562 zze0 = ze0(ji,jj) * EXP( - ze3t * r1_si0 ) ; zzeR = zeR(ji,jj) * EXP( - ze3t * r1_LR ) ! IR ; Red at jk+1 w-level 563 zzeG = zeG(ji,jj) * EXP( - ze3t * r1_LG ) ; zzeB = zeB(ji,jj) * EXP( - ze3t * r1_LB ) ! Green ; Blue - - 564 zzeT = ( zze0 + zzeB + zzeG + zzeR ) * wmask(ji,jj,jk+1) ! Total - - 565 !!st7-9 zzeT = ( zze0 + zzeR + zzeG + zzeB ) * wmask(ji,jj,jk+1) ! Total - - 566 #if defined key_RK3 567 ! ! RK3 : temperature trend at jk t-level 568 pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t 569 #else 570 ! ! MLF : heat content trend due to Qsr flux (qsr_hc) 571 qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) 572 #endif 573 ze0(ji,jj) = zze0 ; zeR(ji,jj) = zzeR ! IR ; Red store at jk+1 w-level 574 zeG(ji,jj) = zzeG ; zeB(ji,jj) = zzeB ! Green ; Blue - - - 575 zeT(ji,jj) = zzeT ! total - - - 576 END_2D 577 !!stbug IF( jk == 1 ) THEN !* sea-ice *! store the 1st level attenuation coeff. 578 !!stbug WHERE( qsr(A2D(0)) /= 0._wp ) ; fraqsr_1lev(A2D(0)) = 1._wp - zeT(A2D(0)) / qsr(A2D(0)) 579 !!stbug ELSEWHERE ; fraqsr_1lev(A2D(0)) = 1._wp 580 !!stbug END WHERE 581 !!stbug ENDIF 582 END DO 583 ! 584 IF( kt == nit000 ) THEN 585 IF(lwp) WRITE(numout,*) 'nk0+1= ', nk0+1, ' qsr IR max = ' , MAXVAL(ze0(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(ze0(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 586 IF(lwp) WRITE(numout,*) ' ', nk0+1, ' qsr R max = ' , MAXVAL(zeR(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeR(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 587 IF(lwp) WRITE(numout,*) ' ', nk0+1, ' qsr G max = ' , MAXVAL(zeG(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeG(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 588 IF(lwp) WRITE(numout,*) ' ', nk0+1, ' qsr B max = ' , MAXVAL(zeB(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeB(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 589 IF(lwp) WRITE(numout,*) ' ', nk0+1, ' qsr T max = ' , MAXVAL(zeT(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeT(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 590 ENDIF 591 ! 592 DO jk = nk0+1, nkR !* down to Red extinction *! (< ~71 meters : RGB , IR removed from calculation) 593 DO_2D( 0, 0, 0, 0 ) 594 ze3t = e3t(ji,jj,jk,Kmm) 595 zzeR = zeR(ji,jj) * EXP( - ze3t * r1_LR ) ! Red at jk+1 w-level 596 zzeG = zeG(ji,jj) * EXP( - ze3t * r1_LG ) ; zzeB = zeB(ji,jj) * EXP( - ze3t * r1_LB ) ! Green ; Blue - - 597 zzeT = ( zzeB + zzeG + zzeR ) * wmask(ji,jj,jk+1) ! Total - - 598 !!st7-11 zzeT = ( zzeR + zzeG + zzeB ) * wmask(ji,jj,jk+1) ! Total - - 599 #if defined key_RK3 600 ! ! RK3 : temperature trend at jk t-level 601 pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t 602 #else 603 ! ! MLF : heat content trend due to Qsr flux (qsr_hc) 604 qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) 605 #endif 606 zeR(ji,jj) = zzeR ! Red store at jk+1 w-level 607 zeG(ji,jj) = zzeG ; zeB(ji,jj) = zzeB ! Green ; Blue - - - 608 zeT(ji,jj) = zzeT ! total - - - 609 END_2D 610 END DO 611 ! 612 IF( kt == nit000 ) THEN 613 IF(lwp) WRITE(numout,*) 'nkR+1= ', nkR+1, ' qsr R max = ' , MAXVAL(zeR(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeR(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 614 IF(lwp) WRITE(numout,*) ' ', nkR+1, ' qsr G max = ' , MAXVAL(zeG(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeG(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 615 IF(lwp) WRITE(numout,*) ' ', nkR+1, ' qsr B max = ' , MAXVAL(zeB(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeB(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 616 IF(lwp) WRITE(numout,*) ' ', nkR+1, ' qsr T max = ' , MAXVAL(zeT(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeT(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 617 ENDIF 618 ! 619 DO jk = nkR+1, nkG !* down to Green extinction *! (< ~350 m : GB , IR+R removed from calculation) 620 DO_2D( 0, 0, 0, 0 ) 621 ze3t = e3t(ji,jj,jk,Kmm) 622 zzeG = zeG(ji,jj) * EXP( - ze3t * r1_LG ) ; zzeB = zeB(ji,jj) * EXP( - ze3t * r1_LB ) ! Green ; Blue at jk+1 w-level 623 zzeT = ( zzeG + zzeB ) * wmask(ji,jj,jk+1) ! Total - - 624 #if defined key_RK3 625 ! ! RK3 : temperature trend at jk t-level 626 pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t 627 #else 628 ! ! MLF : heat content trend due to Qsr flux (qsr_hc) 629 qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) 630 #endif 631 zeG(ji,jj) = zzeG ; zeB(ji,jj) = zzeB ! Green ; Blue store at jk+1 w-level 632 zeT(ji,jj) = zzeT ! total - - - 633 END_2D 634 END DO 635 ! 636 IF( kt == nit000 ) THEN 637 IF(lwp) WRITE(numout,*) 'nkG+1= ', nkG+1, ' qsr G max = ' , MAXVAL(zeG(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeG(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 638 IF(lwp) WRITE(numout,*) ' ', nkG+1, ' qsr B max = ' , MAXVAL(zeB(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeB(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 639 IF(lwp) WRITE(numout,*) ' ', nkG+1, ' qsr T max = ' , MAXVAL(zeT(:,:)*wmask(:,:,jk)), ' W/m2' , MAXVAL(zeT(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)*wmask(:,:,jk)), ' K/s' 640 ENDIF 641 ! 642 DO jk = nkG+1, nkB !* down to Blue extinction *! (< ~1300 m : B , IR+RG removed from calculation) 643 DO_2D( 0, 0, 0, 0 ) 644 ze3t = e3t(ji,jj,jk,Kmm) 645 zzeB = zeB(ji,jj) * EXP( - ze3t * r1_LB ) ! Blue at jk+1 w-level 646 zzeT = ( zzeB ) * wmask(ji,jj,jk+1) ! Total - - 647 #if defined key_RK3 648 ! ! RK3 : temperature trend at jk t-level 649 pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) / ze3t 650 #else 651 ! ! MLF : heat content trend due to Qsr flux (qsr_hc) 652 qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zeT(ji,jj) - zzeT ) 653 #endif 654 zeB(ji,jj) = zzeB ! Blue store at jk+1 w-level 655 zeT(ji,jj) = zzeT ! total - - - 656 END_2D 657 END DO 658 ! 659 IF( kt == nit000 ) THEN 660 IF(lwp) WRITE(numout,*) 'nkB+1= ', nkB+1, ' qsr T max = ' , MAXVAL(zeT), ' W/m2' , MAXVAL(zeT(:,:)*r1_rho0_rcp/e3t(:,:,nk0+1,Kmm)), ' K/s' 661 ENDIF 662 ! 663 END SUBROUTINE qsr_RGB 664 665 666 SUBROUTINE qsr_2BD( Kmm, pts, Krhs ) 667 !!---------------------------------------------------------------------- 668 !! *** ROUTINE qsr_2BD *** 669 !! 670 !! ** Purpose : 2 bands (IR+visible) solar radiation with constant chlorophyll 671 !! 672 !! ** Method : The profile of the solar radiation within the ocean is defined 673 !! through 2 wavebands (rn_si0,rn_si1) a ratio rn_abs for IR absorbtion. 674 !! Considering the 2 wavebands case: 675 !! I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) ) 676 !! The temperature trend associated with the solar radiation penetration 677 !! is given by : zta = 1/e3t dk[ I ] / (rho0*Cp) 678 !! At the bottom, boudary condition for the radiation is no flux : 679 !! all heat which has not been absorbed in the above levels is put 680 !! in the last ocean level. 681 !! The computation is only done down to the level where 682 !! I(k) < 1.e-15 W/m2 (i.e. over the top nk levels) . 683 !! 684 !! ** Action : - update ta with the penetrative solar radiation trend 685 !! - send trend for further diagnostics (l_trdtra=T) 686 !! 687 !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp. 688 !! Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516. 689 !! Morel, A. et Berthon, JF, 1989, Limnol Oceanogr 34(8), 1545-1562 690 !!---------------------------------------------------------------------- 691 INTEGER, INTENT(in ) :: Kmm, Krhs ! ocean time-step and time-level indices 692 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation 693 !! 694 INTEGER :: ji, jj, jk ! dummy loop indices 695 REAL(wp) :: zzatt ! - - 696 REAL(wp) :: zz0 , zz1 , ze3t ! - - 697 REAL(wp), DIMENSION(A2D(0)) :: zatt 698 !!---------------------------------------------------------------------- 699 ! 700 ! !======================! 701 ! !== 2-bands fluxes ==! 702 ! !======================! 703 ! 704 zz0 = rn_abs * r1_rho0_rcp ! surface equi-partition in 2-bands 705 zz1 = ( 1._wp - rn_abs ) * r1_rho0_rcp 706 ! 707 zatt(A2D(0)) = r1_rho0_rcp !* surface value *! 708 ! 709 DO_2D( 0, 0, 0, 0 ) 710 zatt(ji,jj) = ( zz0 * EXP( -gdepw(ji,jj,1,Kmm)*r1_si0 ) + zz1 * EXP( -gdepw(ji,jj,1,Kmm)*r1_si1 ) ) 711 END_2D 712 ! 713 !!st IF(lwp) WRITE(numout,*) 'level = ', 1, ' qsr max = ' , MAXVAL(zatt)*rho0_rcp, ' W/m2', ' qsr min = ' , MINVAL(zatt)*rho0_rcp, ' W/m2' 714 ! 715 DO jk = 1, nk0 !* near surface layers *! (< ~14 meters : IR + visible light ) 716 DO_2D( 0, 0, 0, 0 ) 717 ze3t = e3t(ji,jj,jk,Kmm) ! light attenuation at jk+1 w-level (divided by rho0_rcp) 718 zzatt = ( zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)*r1_si0 ) & 719 & + zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*r1_si1 ) ) * wmask(ji,jj,jk+1) 720 #if defined key_RK3 721 ! ! RK3 : temperature trend at jk t-level 722 pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + qsr(ji,jj) * ( zatt(ji,jj) - zzatt ) / ze3t 723 #else 724 ! ! MLF : heat content trend due to Qsr flux (qsr_hc) 725 qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zatt(ji,jj) - zzatt ) 726 #endif 727 zatt(ji,jj) = zzatt ! save for the next level computation 728 END_2D 729 !!stbug ! !* sea-ice *! store the 1st level attenuation coeff. 730 !!stbug IF( jk == 1 ) fraqsr_1lev(A2D(0)) = 1._wp - zatt(A2D(0)) * rho0_rcp 731 END DO 732 !!st IF(lwp) WRITE(numout,*) 'nk0+1= ', nk0+1, ' qsr max = ' , MAXVAL(zatt*qsr)*rho0_rcp, ' W/m2' , MAXVAL(zatt*qsr/e3t(:,:,nk0+1,Kmm)), ' K/s' 733 ! 734 DO jk = nk0+1, nkV !* deeper layers *! (visible light only) 735 DO_2D( 0, 0, 0, 0 ) 736 ze3t = e3t(ji,jj,jk,Kmm) ! light attenuation at jk+1 w-level (divided by rho0_rcp) 737 zzatt = ( zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*r1_si1 ) ) * wmask(ji,jj,jk+1) 738 #if defined key_RK3 739 ! ! RK3 : temperature trend at jk t-level 740 pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + qsr(ji,jj) * ( zatt(ji,jj) - zzatt ) / ze3t 741 #else 742 ! ! MLF : heat content trend due to Qsr flux (qsr_hc) 743 qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zatt(ji,jj) - zzatt ) 744 #endif 745 zatt(ji,jj) = zzatt ! save for the next level computation 746 END_2D 747 END DO 748 ! 749 !!st IF(lwp) WRITE(numout,*) 'nkV+1= ', nkV+1, ' qsr max = ' , MAXVAL(zatt*qsr)*rho0_rcp, ' W/m2' , MAXVAL(zatt*qsr/e3t(:,:,nkV+1,Kmm)), ' K/s' 750 END SUBROUTINE qsr_2bd 751 752 753 FUNCTION qsr_ext_lev( pL, pfr ) RESULT( klev ) 754 !!---------------------------------------------------------------------- 755 !! *** ROUTINE trc_oce_ext_lev *** 756 !! 757 !! ** Purpose : compute the maximum level of light penetration 758 !! 759 !! ** Method : the function provides the level at which irradiance, I, 760 !! has a negligible effect on temperature. 761 !! T(n+1)-T(n) = ∆t dk[I] / ( rho0 Cp e3t_k ) 762 !! I(k) has a negligible effect on temperature at level k if: 763 !! ∆t I(k) / ( rho0*Cp*e3t_k ) <= 1.e-15 °C 764 !! with I(z) = Qsr*pfr*EXP(-z/L), therefore : 765 !! z >= L * LOG( 1.e-15 * rho0*Cp*e3t_k / ( ∆t*Qsr*pfr ) ) 766 !! with Qsr being the maximum normal surface irradiance at sea 767 !! level (~1000 W/m2). 768 !! # pL is the longest depth of extinction: 769 !! - pL = 23.00 m (2 bands case) 770 !! - pL = 48.24 m (3 bands case: blue waveband & 0.03 mg/m2 for the chlorophyll) 771 !! # pfr is the fraction of solar radiation which penetrates, 772 !! considering Qsr=1000 W/m2 and rn_abs = 0.58: 773 !! - Qsr*pfr0 = Qsr * rn_abs = 580 W/m2 (top absorbtion) 774 !! - Qsr*pfr1 = Qsr * (1-rn_abs) = 420 W/m2 (2 bands case) 775 !! - Qsr*pfr1 = Qsr * (1-rn_abs)/3 = 140 W/m2 (3 bands case & equi-partition) 776 !! 777 !!---------------------------------------------------------------------- 778 INTEGER :: klev ! result: maximum level of light penetration 779 REAL(wp), INTENT(in) :: pL ! depth of extinction 780 REAL(wp), INTENT(in) :: pfr ! frac. solar radiation which penetrates 781 ! 782 INTEGER :: jk ! dummy loop index 783 REAL(wp) :: zcoef ! local scalar 784 REAL(wp) :: zhext ! deepest depth until which light penetrates 785 REAL(wp) :: ze3t , zdw ! max( e3t_k ) and min( w-depth_k+1 ) 786 REAL(wp) :: zprec = 10.e-15_wp ! required precision 787 REAL(wp) :: zQmax= 1000._wp ! maximum normal surface irradiance at sea level (W/m2) 788 !!---------------------------------------------------------------------- 789 ! 790 zQmax = 1000._wp ! maximum normal surface irradiance at sea level (W/m2) 791 ! 792 zcoef = zprec * rho0_rcp / ( rDt * zQmax * pfr) 793 ! 794 klev = jpkm1 ! Level of light extinction 795 DO jk = jpkm1, 1, -1 796 IF( SUM( tmask(:,:,jk) ) > 0 ) THEN ! ocean point at that level 797 zdw = MAXVAL( gdepw_0(:,:,jk+1) * wmask(:,:,jk) ) ! max w-depth at jk+1 level 798 ze3t = MINVAL( e3t_0(:,:,jk ) ) ! minimum e3t at jk level 799 zhext = - pL * LOG( zcoef * ze3t ) ! extinction depth 800 IF( zdw >= zhext ) klev = jk ! last T-level reached by Qsr 801 ELSE ! only land point at level jk 802 klev = jk ! local domain sea-bed level 803 ENDIF 804 END DO 805 ! 806 !!st IF(lwp) WRITE(numout,*) ' level of e3t light extinction = ', klev, ' ref depth = ', gdepw_1d(klev+1), ' m' 807 !!st ! 808 !!st klev = jpkm1 ! Level of light extinction 809 !!st DO jk = jpkm1, 1, -1 810 !!st IF( SUM( tmask(:,:,jk) ) > 0 ) THEN ! ocean point at that level 811 !!st zdw = MAXVAL( gdepw_0(:,:,jk+1) * wmask(:,:,jk) ) ! max w-depth at jk+1 level 812 !!st ze3t = MINVAL( e3t_0(:,:,jk ) ) ! minimum e3t at jk level 813 !!st zhext = - pL * LOG( zcoef * pL ) ! extinction depth 814 !!st IF( zdw >= zhext ) klev = jk ! last T-level reached by Qsr 815 !!st ELSE ! only land point at level jk 816 !!st klev = jk ! local domain sea-bed level 817 !!st ENDIF 818 !!st END DO 819 !!st ! 820 !!st IF(lwp) WRITE(numout,*) ' level of pL light extinction = ', klev, ' ref depth = ', gdepw_1d(klev+1), ' m' 821 ! 822 END FUNCTION qsr_ext_lev 337 823 338 824 … … 355 841 !!---------------------------------------------------------------------- 356 842 INTEGER :: ji, jj, jk ! dummy loop indices 357 INTEGER :: ios, i rgb, ierror, ioptio ! local integer358 REAL(wp) :: z z0, zc0 , zc1, zcoef ! local scalars359 REAL(wp) :: z z1, zc2 , zc3, zchl! - -843 INTEGER :: ios, ierror, ioptio ! local integer 844 REAL(wp) :: zLB, zLG, zLR ! local scalar 845 REAL(wp) :: zVlp, zchl ! - - 360 846 ! 361 847 CHARACTER(len=100) :: cn_dir ! Root directory for location of ssr files … … 368 854 READ ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901) 369 855 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_qsr in reference namelist' ) 370 ! 371 READ ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902 ) 856 READ ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902) 372 857 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist' ) 373 858 IF(lwm) WRITE ( numond, namtra_qsr ) 374 859 ! 375 IF(lwp) THEN ! control print860 IF(lwp) THEN !** control print **! 376 861 WRITE(numout,*) 377 862 WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation' … … 383 868 WRITE(numout,*) ' RGB : Chl data (=1) or cst value (=0) nn_chldta = ', nn_chldta 384 869 WRITE(numout,*) ' RGB & 2 bands: fraction of light (rn_si1) rn_abs = ', rn_abs 385 WRITE(numout,*) ' RGB & 2 bands: shortess depth of extinctionrn_si0 = ', rn_si0386 WRITE(numout,*) ' 2 bands: longest depth of extinctionrn_si1 = ', rn_si1870 WRITE(numout,*) ' RGB & 2 bands: shortess attenuation depth rn_si0 = ', rn_si0 871 WRITE(numout,*) ' 2 bands: longest attenuation depth rn_si1 = ', rn_si1 387 872 WRITE(numout,*) 388 873 ENDIF 389 874 ! 390 ioptio = 0 ! Parameter control875 ioptio = 0 !** Parameter control **! 391 876 IF( ln_qsr_rgb ) ioptio = ioptio + 1 392 877 IF( ln_qsr_2bd ) ioptio = ioptio + 1 … … 401 886 IF( ln_qsr_bio ) nqsr = np_BIO 402 887 ! 403 ! ! Initialisation 404 xsi0r = 1._wp / rn_si0 405 xsi1r = 1._wp / rn_si1 888 ! !** Initialisation **! 889 ! 890 ! !== Infrared attenuation ==! (all schemes) 891 ! !============================! 892 ! 893 r1_si0 = 1._wp / rn_si0 ! inverse of infrared attenuation length 894 ! 895 nk0 = qsr_ext_lev( rn_si0, rn_abs ) ! level of light extinction 896 ! 897 IF(lwp) WRITE(numout,*) ' ==>>> Infrared light attenuation' 898 IF(lwp) WRITE(numout,*) ' level of infrared extinction = ', nk0, ' ref depth = ', gdepw_1d(nk0+1), ' m' 899 IF(lwp) WRITE(numout,*) 406 900 ! 407 901 SELECT CASE( nqsr ) 408 902 ! 409 CASE( np_RGB , np_RGBc ) !== Red-Green-Blue light penetration ==! 410 ! 411 IF(lwp) WRITE(numout,*) ' ==>>> R-G-B light penetration ' 903 CASE( np_RGBc, np_RGB ) !== Red-Green-Blue light attenuation ==! (Chl data or constant) 904 ! !========================================! 905 ! 906 IF( nqsr == np_RGB ) THEN ; zchl = 0.05 ! constant Chl value 907 ELSE ; zchl = 0.03 ! minimum Chl value 908 ENDIF 909 zchl = MAX( 0.03_wp , MIN( zchl , 10._wp) ) ! NB. make sure that chosen value verifies: 0.03 < zchl < 10 910 nc_rgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) ! Convert Chl value to attenuation coefficient look-up table index 412 911 ! 413 912 CALL trc_oce_rgb( rkrgb ) ! tabulated attenuation coef. 414 913 ! 415 nksr = trc_oce_ext_lev( r_si2, 33._wp ) ! level of light extinction 416 ! 417 IF(lwp) WRITE(numout,*) ' level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m' 914 zVlp = ( 1._wp - rn_abs ) / 3._wp ! visible light equi-partition 915 ! 916 ! 1 / length ! attenuation length ! attenuation level 917 r1_LR = rkrgb(3,nc_rgb) ; zLR = 1._wp / r1_LR ; nkR = qsr_ext_lev( zLR, zVlp ) ! Red 918 r1_LG = rkrgb(2,nc_rgb) ; zLG = 1._wp / r1_LG ; nkG = qsr_ext_lev( zLG, zVlp ) ! Green 919 r1_LB = rkrgb(1,nc_rgb) ; zLB = 1._wp / r1_LB ; nkB = qsr_ext_lev( zLB, zVlp ) ! Blue 920 ! 921 nkV = nkB ! maximum level of light penetration 922 ! 923 IF( nqsr == np_RGB ) THEN 924 IF(lwp) WRITE(numout,*) ' ==>>> RGB: light attenuation with a constant Chlorophyll = ', zchl 925 ELSE 926 IF(lwp) WRITE(numout,*) ' ==>>> RGB: light attenuation using Chlorophyll data with min(Chl) = ', zchl 927 ENDIF 928 IF(lwp) WRITE(numout,*) ' level of Red extinction = ', nkR, ' ref depth = ', gdepw_1d(nkR+1), ' m' 929 IF(lwp) WRITE(numout,*) ' level of Green extinction = ', nkG, ' ref depth = ', gdepw_1d(nkG+1), ' m' 930 IF(lwp) WRITE(numout,*) ' level of Blue extinction = ', nkB, ' ref depth = ', gdepw_1d(nkB+1), ' m' 931 IF(lwp) WRITE(numout,*) 418 932 ! 419 933 IF( nqsr == np_RGBc ) THEN ! Chl data : set sf_chl structure … … 423 937 CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_chl structure' ) ; RETURN 424 938 ENDIF 425 ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1) 939 ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1) ) 426 940 IF( sn_chl%ln_tint ) ALLOCATE( sf_chl(1)%fdta(jpi,jpj,1,2) ) 427 941 ! ! fill sf_chl with sn_chl and control print 428 CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init', &942 CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init', & 429 943 & 'Solar penetration function of read chlorophyll', 'namtra_qsr' , no_print ) 430 944 ENDIF 431 IF( nqsr == np_RGB ) THEN ! constant Chl 432 IF(lwp) WRITE(numout,*) ' ==>>> Constant Chlorophyll concentration = 0.05' 433 ENDIF 434 ! 435 CASE( np_2BD ) !== 2 bands light penetration ==! 436 ! 437 IF(lwp) WRITE(numout,*) ' ==>>> 2 bands light penetration' 438 ! 439 nksr = trc_oce_ext_lev( rn_si1, 100._wp ) ! level of light extinction 440 IF(lwp) WRITE(numout,*) ' level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m' 945 ! 946 CASE( np_2BD ) !== 2 bands light attenuation (IR+ visible light) ==! 947 ! 948 ! 949 r1_si1 = 1._wp / rn_si1 ! inverse of visible light attenuation 950 zVlp = ( 1._wp - rn_abs ) ! visible light partition 951 nkV = qsr_ext_lev( rn_si1, zVlp ) ! level of visible light extinction 952 ! 953 IF(lwp) WRITE(numout,*) ' ==>>> 2 bands attenuation (Infrared + Visible light) ' 954 IF(lwp) WRITE(numout,*) ' level of visible light extinction = ', nkV, ' ref depth = ', gdepw_1d(nkV+1), ' m' 955 IF(lwp) WRITE(numout,*) 441 956 ! 442 957 CASE( np_BIO ) !== BIO light penetration ==! … … 447 962 CALL trc_oce_rgb( rkrgb ) ! tabulated attenuation coef. 448 963 ! 449 nk sr = trc_oce_ext_lev( r_si2, 33._wp ) !level of light extinction450 ! 451 IF(lwp) WRITE(numout,*) ' level of light extinction = ', nk sr, ' ref depth = ', gdepw_1d(nksr+1), ' m'964 nkV = trc_oce_ext_lev( r_si2, 33._wp ) ! maximum level of light extinction 965 ! 966 IF(lwp) WRITE(numout,*) ' level of light extinction = ', nkV, ' ref depth = ', gdepw_1d(nkV+1), ' m' 452 967 ! 453 968 END SELECT 454 969 ! 455 qsr_hc(:,:,:) = 0._wp ! now qsr heat content set to zero where it will not be computed 456 ! 457 ! 1st ocean level attenuation coefficient (used in sbcssm) 970 nksr = nkV ! name of max level of light extinction used in traatf(_qco).F90 971 ! 972 #if ! defined key_RK3 973 qsr_hc(:,:,:) = 0._wp ! MLF : now qsr heat content set to zero where it will not be computed 974 #endif 975 ! 976 ! ! Sea-ice : 1st ocean level attenuation coefficient (used in sbcssm) 458 977 IF( iom_varid( numror, 'fraqsr_1lev', ldstop = .FALSE. ) > 0 ) THEN 459 978 CALL iom_get( numror, jpdom_auto, 'fraqsr_1lev' , fraqsr_1lev )
Note: See TracChangeset
for help on using the changeset viewer.
