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