Changeset 6140 for trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
- Timestamp:
- 2015-12-21T12:35:23+01:00 (8 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90 (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
r5836 r6140 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 DO jj = 1, jpj ! Separation in R-G-B depending of the surface Chl 192 DO ji = 1, jpi 193 zchl = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) 194 irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) 195 zekb(ji,jj) = rkrgb(1,irgb) 196 zekg(ji,jj) = rkrgb(2,irgb) 197 zekr(ji,jj) = rkrgb(3,irgb) 198 END DO 199 END DO 200 ELSE ! Variable ocean volume but constant chrlorophyll 201 zchl = 0.05 ! constant chlorophyll 202 irgb = NINT( 41 + 20.*LOG10( zchl ) + 1.e-15 ) 203 zekb(:,:) = rkrgb(1,irgb) ! Separation in R-G-B depending of the chlorophyll 204 zekg(:,:) = rkrgb(2,irgb) 205 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 206 255 ENDIF 207 ! 208 zcoef = ( 1. - rn_abs ) / 3.e0 ! equi-partition in R-G-B 209 ze0(:,:,1) = rn_abs * qsr(:,:) 210 ze1(:,:,1) = zcoef * qsr(:,:) 211 ze2(:,:,1) = zcoef * qsr(:,:) 212 ze3(:,:,1) = zcoef * qsr(:,:) 213 zea(:,:,1) = qsr(:,:) 214 ! 215 DO jk = 2, nksr+1 216 DO jj = 1, jpj 217 DO ji = 1, jpi 218 zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r ) 219 zc1 = ze1(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekb(ji,jj) ) 220 zc2 = ze2(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekg(ji,jj) ) 221 zc3 = ze3(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekr(ji,jj) ) 222 ze0(ji,jj,jk) = zc0 223 ze1(ji,jj,jk) = zc1 224 ze2(ji,jj,jk) = zc2 225 ze3(ji,jj,jk) = zc3 226 zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,jk) 227 END DO 228 END DO 229 END DO 230 ! clem: store attenuation coefficient of the first ocean level 231 IF ( ln_qsr_ice ) THEN 232 DO jj = 1, jpj 233 DO ji = 1, jpi 234 zzc0 = rn_abs * EXP( - fse3t(ji,jj,1) * xsi0r ) 235 zzc1 = zcoef * EXP( - fse3t(ji,jj,1) * zekb(ji,jj) ) 236 zzc2 = zcoef * EXP( - fse3t(ji,jj,1) * zekg(ji,jj) ) 237 zzc3 = zcoef * EXP( - fse3t(ji,jj,1) * zekr(ji,jj) ) 238 fraqsr_1lev(ji,jj) = 1.0 - ( zzc0 + zzc1 + zzc2 + zzc3 ) * tmask(ji,jj,2) 239 END DO 240 END DO 241 ENDIF 242 ! 243 DO jk = 1, nksr ! compute and add qsr trend to ta 244 qsr_hc(:,:,jk) = r1_rau0_rcp * ( zea(:,:,jk) - zea(:,:,jk+1) ) 245 END DO 246 zea(:,:,nksr+1:jpk) = 0.e0 ! below 400m set to zero 247 CALL iom_put( 'qsr3d', zea ) ! Shortwave Radiation 3D distribution 248 ! 249 ELSE !* Constant Chlorophyll 250 DO jk = 1, nksr 251 qsr_hc(:,:,jk) = etot3(:,:,jk) * qsr(:,:) 252 END DO 253 ! clem: store attenuation coefficient of the first ocean level 254 IF ( ln_qsr_ice ) THEN 255 fraqsr_1lev(:,:) = etot3(:,:,1) / r1_rau0_rcp 256 ENDIF 257 ENDIF 258 259 ENDIF 260 ! ! ------------------------- ! 261 IF( ln_qsr_2bd ) THEN ! 2 band light penetration ! 262 ! ! ------------------------- ! 263 ! 264 IF( lk_vvl ) THEN !* variable volume 265 zz0 = rn_abs * r1_rau0_rcp 266 zz1 = ( 1. - rn_abs ) * r1_rau0_rcp 267 DO jk = 1, nksr ! solar heat absorbed at T-point in the top 400m 268 DO jj = 1, jpj 269 DO ji = 1, jpi 270 zc0 = zz0 * EXP( -fsdepw(ji,jj,jk )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk )*xsi1r ) 271 zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r ) 272 qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) ) 273 END DO 274 END DO 275 END DO 276 ! clem: store attenuation coefficient of the first ocean level 277 IF ( ln_qsr_ice ) THEN 278 DO jj = 1, jpj 279 DO ji = 1, jpi 280 zc0 = zz0 * EXP( -fsdepw(ji,jj,1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,1)*xsi1r ) 281 zc1 = zz0 * EXP( -fsdepw(ji,jj,2)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,2)*xsi1r ) 282 fraqsr_1lev(ji,jj) = ( zc0*tmask(ji,jj,1) - zc1*tmask(ji,jj,2) ) / r1_rau0_rcp 283 END DO 284 END DO 285 ENDIF 286 ELSE !* constant volume: coef. computed one for all 287 DO jk = 1, nksr 288 DO jj = 2, jpjm1 289 DO ji = fs_2, fs_jpim1 ! vector opt. 290 ! (ISF) no light penetration below the ice shelves 291 qsr_hc(ji,jj,jk) = etot3(ji,jj,jk) * qsr(ji,jj) * tmask(ji,jj,1) 292 END DO 293 END DO 294 END DO 295 ! clem: store attenuation coefficient of the first ocean level 296 IF ( ln_qsr_ice ) THEN 297 fraqsr_1lev(:,:) = etot3(:,:,1) / r1_rau0_rcp 298 ENDIF 299 ! 300 ENDIF 301 ! 302 ENDIF 303 ! 304 ! Add to the general trend 305 DO jk = 1, nksr 306 DO jj = 2, jpjm1 307 DO ji = fs_2, fs_jpim1 ! vector opt. 308 z1_e3t = zfact / fse3t(ji,jj,jk) 309 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 310 END DO 311 END DO 312 END DO 313 ! 314 ENDIF 315 ! 316 IF( lrst_oce ) THEN ! Write in the ocean restart file 317 ! ******************************* 318 IF(lwp) WRITE(numout,*) 319 IF(lwp) WRITE(numout,*) 'qsr tracer content forcing field written in ocean restart file ', & 320 & 'at it= ', kt,' date= ', ndastp 321 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 322 275 CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b' , qsr_hc ) 323 CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev ) ! default definition in sbcssm 324 ! 325 ENDIF 326 276 CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev ) 277 ENDIF 278 ! 327 279 IF( l_trdtra ) THEN ! qsr tracers trends saved for diagnostics 328 280 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 329 281 CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt ) 330 CALL wrk_dealloc( jpi, jpj, jpk,ztrdt )282 CALL wrk_dealloc( jpi,jpj,jpk, ztrdt ) 331 283 ENDIF 332 284 ! ! print mean trends (used for debugging) 333 285 IF(ln_ctl) CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' qsr - Ta: ', mask1=tmask, clinfo3='tra-ta' ) 334 !335 CALL wrk_dealloc( jpi, jpj, zekb, zekg, zekr )336 CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )337 286 ! 338 287 IF( nn_timing == 1 ) CALL timing_stop('tra_qsr') … … 358 307 !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp. 359 308 !!---------------------------------------------------------------------- 360 ! 361 INTEGER :: ji, jj, jk ! dummy loop indices 362 INTEGER :: irgb, ierror, ioptio, nqsr ! local integer 363 INTEGER :: ios ! Local integer output status for namelist read 364 REAL(wp) :: zz0, zc0 , zc1, zcoef ! local scalars 365 REAL(wp) :: zz1, zc2 , zc3, zchl ! - - 366 REAL(wp), POINTER, DIMENSION(:,: ) :: zekb, zekg, zekr 367 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 ! - - 368 313 ! 369 314 CHARACTER(len=100) :: cn_dir ! Root directory for location of ssr files 370 315 TYPE(FLD_N) :: sn_chl ! informations about the chlorofyl field to be read 371 316 !! 372 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, & 373 318 & nn_chldta, rn_abs, rn_si0, rn_si1 374 319 !!---------------------------------------------------------------------- 375 376 ! 377 IF( nn_timing == 1 ) CALL timing_start('tra_qsr_init') 378 ! 379 CALL wrk_alloc( jpi, jpj, zekb, zekg, zekr ) 380 CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 381 ! 382 383 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 384 324 READ ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901) 385 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_qsr in reference namelist', lwp )386 387 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 388 328 READ ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902 ) 389 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 ) 390 330 IF(lwm) WRITE ( numond, namtra_qsr ) 391 331 ! … … 395 335 WRITE(numout,*) '~~~~~~~~~~~~' 396 336 WRITE(numout,*) ' Namelist namtra_qsr : set the parameter of penetration' 397 WRITE(numout,*) ' Light penetration (T) or not (F) ln_traqsr = ', ln_traqsr 398 WRITE(numout,*) ' RGB (Red-Green-Blue) light penetration ln_qsr_rgb = ', ln_qsr_rgb 399 WRITE(numout,*) ' 2 band light penetration ln_qsr_2bd = ', ln_qsr_2bd 400 WRITE(numout,*) ' bio-model light penetration ln_qsr_bio = ', ln_qsr_bio 401 WRITE(numout,*) ' light penetration for ice-model LIM3 ln_qsr_ice = ', ln_qsr_ice 402 WRITE(numout,*) ' RGB : Chl data (=1) or cst value (=0) nn_chldta = ', nn_chldta 403 WRITE(numout,*) ' RGB & 2 bands: fraction of light (rn_si1) rn_abs = ', rn_abs 404 WRITE(numout,*) ' RGB & 2 bands: shortess depth of extinction rn_si0 = ', rn_si0 405 WRITE(numout,*) ' 2 bands: longest depth of extinction rn_si1 = ', rn_si1 406 ENDIF 407 408 IF( ln_traqsr ) THEN ! control consistency 409 ! 410 IF( .NOT.lk_qsr_bio .AND. ln_qsr_bio ) THEN 411 CALL ctl_warn( 'No bio model : force ln_qsr_bio = FALSE ' ) 412 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' ) 413 388 ENDIF 414 ! 415 ioptio = 0 ! Parameter control 416 IF( ln_qsr_rgb ) ioptio = ioptio + 1 417 IF( ln_qsr_2bd ) ioptio = ioptio + 1 418 IF( ln_qsr_bio ) ioptio = ioptio + 1 419 ! 420 IF( ioptio /= 1 ) & 421 CALL ctl_stop( ' Choose ONE type of light penetration in namelist namtra_qsr', & 422 & ' 2 bands, 3 RGB bands or bio-model light penetration' ) 423 ! 424 IF( ln_qsr_rgb .AND. nn_chldta == 0 ) nqsr = 1 425 IF( ln_qsr_rgb .AND. nn_chldta == 1 ) nqsr = 2 426 IF( ln_qsr_2bd ) nqsr = 3 427 IF( ln_qsr_bio ) nqsr = 4 428 ! 429 IF(lwp) THEN ! Print the choice 430 WRITE(numout,*) 431 IF( nqsr == 1 ) WRITE(numout,*) ' R-G-B light penetration - Constant Chlorophyll' 432 IF( nqsr == 2 ) WRITE(numout,*) ' R-G-B light penetration - Chl data ' 433 IF( nqsr == 3 ) WRITE(numout,*) ' 2 bands light penetration' 434 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' 435 391 ENDIF 436 392 ! 437 ENDIF 438 ! ! ===================================== ! 439 IF( ln_traqsr ) THEN ! Initialisation of Light Penetration ! 440 ! ! ===================================== ! 441 ! 442 xsi0r = 1.e0 / rn_si0 443 xsi1r = 1.e0 / rn_si1 444 ! ! ---------------------------------- ! 445 IF( ln_qsr_rgb ) THEN ! Red-Green-Blue light penetration ! 446 ! ! ---------------------------------- ! 447 ! 448 CALL trc_oce_rgb( rkrgb ) !* tabulated attenuation coef. 449 ! 450 ! !* level of light extinction 451 IF( ln_sco ) THEN ; nksr = jpkm1 452 ELSE ; nksr = trc_oce_ext_lev( r_si2, 0.33e2 ) 453 ENDIF 454 455 IF(lwp) WRITE(numout,*) ' level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m' 456 ! 457 IF( nn_chldta == 1 ) THEN !* Chl data : set sf_chl structure 458 IF(lwp) WRITE(numout,*) 459 IF(lwp) WRITE(numout,*) ' Chlorophyll read in a file' 460 ALLOCATE( sf_chl(1), STAT=ierror ) 461 IF( ierror > 0 ) THEN 462 CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_chl structure' ) ; RETURN 463 ENDIF 464 ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1) ) 465 IF( sn_chl%ln_tint )ALLOCATE( sf_chl(1)%fdta(jpi,jpj,1,2) ) 466 ! ! fill sf_chl with sn_chl and control print 467 CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init', & 468 & 'Solar penetration function of read chlorophyll', 'namtra_qsr' ) 469 ! 470 ELSE !* constant Chl : compute once for all the distribution of light (etot3) 471 IF(lwp) WRITE(numout,*) 472 IF(lwp) WRITE(numout,*) ' Constant Chlorophyll concentration = 0.05' 473 IF( lk_vvl ) THEN ! variable volume 474 IF(lwp) WRITE(numout,*) ' key_vvl: light distribution will be computed at each time step' 475 ELSE ! constant volume: computes one for all 476 IF(lwp) WRITE(numout,*) ' fixed volume: light distribution computed one for all' 477 ! 478 zchl = 0.05 ! constant chlorophyll 479 irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) 480 zekb(:,:) = rkrgb(1,irgb) ! Separation in R-G-B depending of the chlorophyll 481 zekg(:,:) = rkrgb(2,irgb) 482 zekr(:,:) = rkrgb(3,irgb) 483 ! 484 zcoef = ( 1. - rn_abs ) / 3.e0 ! equi-partition in R-G-B 485 ze0(:,:,1) = rn_abs 486 ze1(:,:,1) = zcoef 487 ze2(:,:,1) = zcoef 488 ze3(:,:,1) = zcoef 489 zea(:,:,1) = tmask(:,:,1) ! = ( ze0+ze1+z2+ze3 ) * tmask 490 491 DO jk = 2, nksr+1 492 DO jj = 1, jpj 493 DO ji = 1, jpi 494 zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * xsi0r ) 495 zc1 = ze1(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekb(ji,jj) ) 496 zc2 = ze2(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekg(ji,jj) ) 497 zc3 = ze3(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekr(ji,jj) ) 498 ze0(ji,jj,jk) = zc0 499 ze1(ji,jj,jk) = zc1 500 ze2(ji,jj,jk) = zc2 501 ze3(ji,jj,jk) = zc3 502 zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,jk) 503 END DO 504 END DO 505 END DO 506 ! 507 DO jk = 1, nksr 508 ! (ISF) no light penetration below the ice shelves 509 etot3(:,:,jk) = r1_rau0_rcp * ( zea(:,:,jk) - zea(:,:,jk+1) ) * tmask(:,:,1) 510 END DO 511 etot3(:,:,nksr+1:jpk) = 0.e0 ! below 400m set to zero 512 ENDIF 513 ENDIF 514 ! 515 ENDIF 516 ! ! ---------------------------------- ! 517 IF( ln_qsr_2bd ) THEN ! 2 bands light penetration ! 518 ! ! ---------------------------------- ! 519 ! 520 ! ! level of light extinction 521 nksr = trc_oce_ext_lev( rn_si1, 1.e2 ) 522 IF(lwp) THEN 523 WRITE(numout,*) 524 IF(lwp) WRITE(numout,*) ' level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m' 525 ENDIF 526 ! 527 IF( lk_vvl ) THEN ! variable volume 528 IF(lwp) WRITE(numout,*) ' key_vvl: light distribution will be computed at each time step' 529 ELSE ! constant volume: computes one for all 530 zz0 = rn_abs * r1_rau0_rcp 531 zz1 = ( 1. - rn_abs ) * r1_rau0_rcp 532 DO jk = 1, nksr !* solar heat absorbed at T-point computed once for all 533 DO jj = 1, jpj ! top 400 meters 534 DO ji = 1, jpi 535 zc0 = zz0 * EXP( -fsdepw(ji,jj,jk )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk )*xsi1r ) 536 zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r ) 537 etot3(ji,jj,jk) = ( zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1) ) * tmask(ji,jj,1) 538 END DO 539 END DO 540 END DO 541 etot3(:,:,nksr+1:jpk) = 0.e0 ! below 400m set to zero 542 ! 543 ENDIF 544 ENDIF 545 ! ! ===================================== ! 546 ELSE ! No light penetration ! 547 ! ! ===================================== ! 548 IF(lwp) THEN 549 WRITE(numout,*) 550 WRITE(numout,*) 'tra_qsr_init : NO solar flux penetration' 551 WRITE(numout,*) '~~~~~~~~~~~~' 552 ENDIF 553 ENDIF 554 ! 555 ! 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) 556 410 IF( iom_varid( numror, 'fraqsr_1lev', ldstop = .FALSE. ) > 0 ) THEN 557 411 CALL iom_get( numror, jpdom_autoglo, 'fraqsr_1lev' , fraqsr_1lev ) 558 412 ELSE 559 fraqsr_1lev(:,:) = 1._wp ! default definition 560 ENDIF 561 ! 562 CALL wrk_dealloc( jpi, jpj, zekb, zekg, zekr ) 563 CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 564 ! 565 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') 566 417 ! 567 418 END SUBROUTINE tra_qsr_init
Note: See TracChangeset
for help on using the changeset viewer.
