- Timestamp:
- 2016-01-08T10:35:19+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
r4624 r6225 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 23 USE trdmod_oce ! ocean variables trends 24 USE trdtra ! ocean active tracers trends 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 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 25 USE trd_oce ! trends: ocean variables 26 USE trdtra ! trends manager: tracers 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 34 USE sbc_ice, ONLY : lk_lim335 37 36 38 IMPLICIT NONE … … 38 40 39 41 PUBLIC tra_qsr ! routine called by step.F90 (ln_traqsr=T) 40 PUBLIC tra_qsr_init ! routine called by opa.F9042 PUBLIC tra_qsr_init ! routine called by nemogcm.F90 41 43 42 44 ! !!* Namelist namtra_qsr: penetrative solar radiation … … 50 52 REAL(wp), PUBLIC :: rn_si0 !: very near surface depth of extinction (RGB & 2 bands) 51 53 REAL(wp), PUBLIC :: rn_si1 !: deepest depth of extinction (water type I) (2 bands) 52 53 ! Module variables 54 REAL(wp) :: xsi0r !: inverse of rn_si0 55 REAL(wp) :: xsi1r !: inverse of rn_si1 54 ! 55 INTEGER , PUBLIC :: nksr !: levels below which the light cannot penetrate (depth larger than 391 m) 56 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 56 67 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_chl ! structure of input Chl (file informations, fields read) 57 INTEGER, PUBLIC :: nksr ! levels below which the light cannot penetrate ( depth larger than 391 m)58 REAL(wp), DIMENSION(3,61) :: rkrgb !: tabulated attenuation coefficients for RGB absorption59 68 60 69 !! * Substitutions 61 # include "domzgr_substitute.h90"62 70 # include "vectopt_loop_substitute.h90" 63 71 !!---------------------------------------------------------------------- … … 73 81 !! 74 82 !! ** Purpose : Compute the temperature trend due to the solar radiation 75 !! penetration and add it to the general temperature trend.83 !! penetration and add it to the general temperature trend. 76 84 !! 77 85 !! ** Method : The profile of the solar radiation within the ocean is defined … … 84 92 !! all heat which has not been absorbed in the above levels is put 85 93 !! in the last ocean level. 86 !! In z-coordinate case, the computation is only done down to the 87 !! level where I(k) < 1.e-15 W/m2. In addition, the coefficients 88 !! used for the computation are calculated one for once as they 89 !! depends on k only. 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) . 90 96 !! 91 97 !! ** Action : - update ta with the penetrative solar radiation trend 92 !! - s ave the trend in ttrd ('key_trdtra')98 !! - send trend for further diagnostics (l_trdtra=T) 93 99 !! 94 100 !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp. 95 101 !! Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516. 96 102 !!---------------------------------------------------------------------- 97 !98 103 INTEGER, INTENT(in) :: kt ! ocean time-step 99 104 ! 100 INTEGER :: ji, jj, jk ! dummy loop indices101 INTEGER :: irgb ! local integers102 REAL(wp) :: zchl, zcoef, z fact! local scalars103 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 ! - - 104 109 REAL(wp) :: zzc0, zzc1, zzc2, zzc3 ! - - 105 REAL(wp) :: zz0 , zz1, z1_e3t! - -106 REAL(wp), POINTER, DIMENSION(:,: ):: zekb, zekg, zekr110 REAL(wp) :: zz0 , zz1 ! - - 111 REAL(wp), POINTER, DIMENSION(:,:) :: zekb, zekg, zekr 107 112 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt 113 REAL(wp), POINTER, DIMENSION(:,:,:) :: zetot 108 114 !!---------------------------------------------------------------------- 109 115 ! 110 116 IF( nn_timing == 1 ) CALL timing_start('tra_qsr') 111 !112 CALL wrk_alloc( jpi, jpj, zekb, zekg, zekr )113 CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )114 117 ! 115 118 IF( kt == nit000 ) THEN … … 117 120 IF(lwp) WRITE(numout,*) 'tra_qsr : penetration of the surface solar radiation' 118 121 IF(lwp) WRITE(numout,*) '~~~~~~~' 119 IF( .NOT.ln_traqsr ) RETURN 120 ENDIF 121 122 IF( l_trdtra ) THEN ! Save ta and sa trends 123 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 ) 124 126 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 125 127 ENDIF 126 127 ! Set before qsr tracer content field128 ! ***********************************129 IF( kt == nit000 ) THEN ! Set the forcing field at nit000 - 1130 ! ! -----------------------------------131 IF( ln_rstart .AND. & ! Restart: read in restart file 132 & iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN133 IF(lwp) WRITE(numout,*) ' nit000-1 qsr tracer content forcing field re d in the restart file'134 z fact = 0.5e0128 ! 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 135 137 CALL iom_get( numror, jpdom_autoglo, 'qsr_hc_b', qsr_hc_b ) ! before heat content trend due to Qsr flux 136 138 ELSE ! No restart or restart not found: Euler forward time stepping 137 z fact = 1.e0138 qsr_hc_b(:,:,:) = 0. e0139 z1_2 = 1._wp 140 qsr_hc_b(:,:,:) = 0._wp 139 141 ENDIF 140 ELSE ! Swap of forcing field 141 ! ! --------------------- 142 zfact = 0.5e0 142 ELSE !== Swap of qsr heat content ==! 143 z1_2 = 0.5_wp 143 144 qsr_hc_b(:,:,:) = qsr_hc(:,:,:) 144 145 ENDIF 145 ! Compute now qsr tracer content field 146 ! ************************************ 147 148 ! ! ============================================== ! 149 IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN ! bio-model fluxes : all vertical coordinates ! 150 ! ! ============================================== ! 151 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 152 154 qsr_hc(:,:,jk) = r1_rau0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) ) 153 155 END DO 154 ! Add to the general trend 155 DO jk = 1, jpkm1 156 DO jj = 2, jpjm1 157 DO ji = fs_2, fs_jpim1 ! vector opt. 158 z1_e3t = zfact / fse3t(ji,jj,jk) 159 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) * z1_e3t 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) 160 171 END DO 161 172 END DO 162 END DO 163 CALL iom_put( 'qsr3d', etot3 ) ! Shortwave Radiation 3D distribution 164 ! clem: store attenuation coefficient of the first ocean level 165 IF ( lk_lim3 .AND. ln_qsr_ice ) THEN 166 DO jj = 1, jpj 167 DO ji = 1, jpi 168 IF ( qsr(ji,jj) /= 0._wp ) THEN 169 oatte(ji,jj) = ( qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) ) ) 170 iatte(ji,jj) = oatte(ji,jj) 171 ENDIF 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) 172 182 END DO 173 183 END DO 174 184 ENDIF 175 ! ! ============================================== ! 176 ELSE ! Ocean alone : 177 ! ! ============================================== ! 178 ! 179 ! ! ------------------------- ! 180 IF( ln_qsr_rgb) THEN ! R-G-B light penetration ! 181 ! ! ------------------------- ! 182 ! Set chlorophyl concentration 183 IF( nn_chldta == 1 .OR. lk_vvl ) THEN !* Variable Chlorophyll or ocean volume 184 ! 185 IF( nn_chldta == 1 ) THEN !* Variable Chlorophyll 186 ! 187 CALL fld_read( kt, 1, sf_chl ) ! Read Chl data and provides it at the current time step 188 ! 189 !CDIR COLLAPSE 190 !CDIR NOVERRCHK 191 DO jj = 1, jpj ! Separation in R-G-B depending of the surface Chl 192 !CDIR NOVERRCHK 193 DO ji = 1, jpi 194 zchl = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) 195 irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) 196 zekb(ji,jj) = rkrgb(1,irgb) 197 zekg(ji,jj) = rkrgb(2,irgb) 198 zekr(ji,jj) = rkrgb(3,irgb) 199 END DO 200 END DO 201 ELSE ! Variable ocean volume but constant chrlorophyll 202 zchl = 0.05 ! constant chlorophyll 203 irgb = NINT( 41 + 20.*LOG10( zchl ) + 1.e-15 ) 204 zekb(:,:) = rkrgb(1,irgb) ! Separation in R-G-B depending of the chlorophyll 205 zekg(:,:) = rkrgb(2,irgb) 206 zekr(:,:) = rkrgb(3,irgb) 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 207 255 ENDIF 208 ! 209 zcoef = ( 1. - rn_abs ) / 3.e0 ! equi-partition in R-G-B 210 ze0(:,:,1) = rn_abs * qsr(:,:) 211 ze1(:,:,1) = zcoef * qsr(:,:) 212 ze2(:,:,1) = zcoef * qsr(:,:) 213 ze3(:,:,1) = zcoef * qsr(:,:) 214 zea(:,:,1) = qsr(:,:) 215 ! 216 DO jk = 2, nksr+1 217 !CDIR NOVERRCHK 218 DO jj = 1, jpj 219 !CDIR NOVERRCHK 220 DO ji = 1, jpi 221 zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r ) 222 zc1 = ze1(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekb(ji,jj) ) 223 zc2 = ze2(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekg(ji,jj) ) 224 zc3 = ze3(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekr(ji,jj) ) 225 ze0(ji,jj,jk) = zc0 226 ze1(ji,jj,jk) = zc1 227 ze2(ji,jj,jk) = zc2 228 ze3(ji,jj,jk) = zc3 229 zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,jk) 230 END DO 231 END DO 232 END DO 233 ! clem: store attenuation coefficient of the first ocean level 234 IF ( lk_lim3 .AND. ln_qsr_ice ) THEN 235 DO jj = 1, jpj 236 DO ji = 1, jpi 237 zzc0 = rn_abs * EXP( - fse3t(ji,jj,1) * xsi0r ) 238 zzc1 = zcoef * EXP( - fse3t(ji,jj,1) * zekb(ji,jj) ) 239 zzc2 = zcoef * EXP( - fse3t(ji,jj,1) * zekg(ji,jj) ) 240 zzc3 = zcoef * EXP( - fse3t(ji,jj,1) * zekr(ji,jj) ) 241 oatte(ji,jj) = 1.0 - ( zzc0 + zzc1 + zzc2 + zzc3 ) * tmask(ji,jj,2) 242 iatte(ji,jj) = 1.0 - ( zzc0 + zzc1 + zcoef + zcoef ) * tmask(ji,jj,2) 243 END DO 244 END DO 245 ENDIF 246 ! 247 DO jk = 1, nksr ! compute and add qsr trend to ta 248 qsr_hc(:,:,jk) = r1_rau0_rcp * ( zea(:,:,jk) - zea(:,:,jk+1) ) 249 END DO 250 zea(:,:,nksr+1:jpk) = 0.e0 ! below 400m set to zero 251 CALL iom_put( 'qsr3d', zea ) ! Shortwave Radiation 3D distribution 252 ! 253 ELSE !* Constant Chlorophyll 254 DO jk = 1, nksr 255 qsr_hc(:,:,jk) = etot3(:,:,jk) * qsr(:,:) 256 END DO 257 ! clem: store attenuation coefficient of the first ocean level 258 IF ( lk_lim3 .AND. ln_qsr_ice ) THEN 259 oatte(:,:) = etot3(:,:,1) / r1_rau0_rcp 260 iatte(:,:) = oatte(:,:) 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 ( lk_lim3 .AND. 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 oatte(ji,jj) = ( zc0*tmask(ji,jj,1) - zc1*tmask(ji,jj,2) ) / r1_rau0_rcp 288 iatte(ji,jj) = oatte(ji,jj) 289 END DO 290 END DO 291 ENDIF 292 ELSE !* constant volume: coef. computed one for all 293 DO jk = 1, nksr 294 DO jj = 2, jpjm1 295 DO ji = fs_2, fs_jpim1 ! vector opt. 296 qsr_hc(ji,jj,jk) = etot3(ji,jj,jk) * qsr(ji,jj) 297 END DO 298 END DO 299 END DO 300 ! clem: store attenuation coefficient of the first ocean level 301 IF ( lk_lim3 .AND. ln_qsr_ice ) THEN 302 oatte(:,:) = etot3(:,:,1) / r1_rau0_rcp 303 iatte(:,:) = oatte(:,:) 304 ENDIF 305 ! 306 ENDIF 307 ! 308 ENDIF 309 ! 310 ! Add to the general trend 311 DO jk = 1, nksr 312 DO jj = 2, jpjm1 313 DO ji = fs_2, fs_jpim1 ! vector opt. 314 z1_e3t = zfact / fse3t(ji,jj,jk) 315 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) + ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) * z1_e3t 316 END DO 317 END DO 318 END DO 319 ! 320 ENDIF 321 ! 322 IF( lrst_oce ) THEN ! Write in the ocean restart file 323 ! ******************************* 324 IF(lwp) WRITE(numout,*) 325 IF(lwp) WRITE(numout,*) 'qsr tracer content forcing field written in ocean restart file ', & 326 & 'at it= ', kt,' date= ', ndastp 327 IF(lwp) WRITE(numout,*) '~~~~' 328 CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b', qsr_hc ) 329 ! 330 ENDIF 331 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 275 CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b' , qsr_hc ) 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 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ trd_qsr, ztrdt )335 CALL wrk_dealloc( jpi, jpj, jpk,ztrdt )281 CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, 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 ! clem init for oatte and iatte 385 IF( .NOT. ln_rstart ) THEN 386 oatte(:,:) = 1._wp 387 iatte(:,:) = 1._wp 388 ENDIF 389 ! 390 CALL wrk_alloc( jpi, jpj, zekb, zekg, zekr ) 391 CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 392 ! 393 394 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 395 324 READ ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901) 396 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_qsr in reference namelist', lwp )397 398 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 399 328 READ ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902 ) 400 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 ) 401 330 IF(lwm) WRITE ( numond, namtra_qsr ) 402 331 ! … … 406 335 WRITE(numout,*) '~~~~~~~~~~~~' 407 336 WRITE(numout,*) ' Namelist namtra_qsr : set the parameter of penetration' 408 WRITE(numout,*) ' Light penetration (T) or not (F) ln_traqsr = ', ln_traqsr 409 WRITE(numout,*) ' RGB (Red-Green-Blue) light penetration ln_qsr_rgb = ', ln_qsr_rgb 410 WRITE(numout,*) ' 2 band light penetration ln_qsr_2bd = ', ln_qsr_2bd 411 WRITE(numout,*) ' bio-model light penetration ln_qsr_bio = ', ln_qsr_bio 412 WRITE(numout,*) ' light penetration for ice-model LIM3 ln_qsr_ice = ', ln_qsr_ice 413 WRITE(numout,*) ' RGB : Chl data (=1) or cst value (=0) nn_chldta = ', nn_chldta 414 WRITE(numout,*) ' RGB & 2 bands: fraction of light (rn_si1) rn_abs = ', rn_abs 415 WRITE(numout,*) ' RGB & 2 bands: shortess depth of extinction rn_si0 = ', rn_si0 416 WRITE(numout,*) ' 2 bands: longest depth of extinction rn_si1 = ', rn_si1 417 WRITE(numout,*) ' light penetration for ice-model LIM3 ln_qsr_ice = ', ln_qsr_ice 418 ENDIF 419 420 IF( ln_traqsr ) THEN ! control consistency 421 ! 422 IF( .NOT.lk_qsr_bio .AND. ln_qsr_bio ) THEN 423 CALL ctl_warn( 'No bio model : force ln_qsr_bio = FALSE ' ) 424 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' ) 425 388 ENDIF 426 ! 427 ioptio = 0 ! Parameter control 428 IF( ln_qsr_rgb ) ioptio = ioptio + 1 429 IF( ln_qsr_2bd ) ioptio = ioptio + 1 430 IF( ln_qsr_bio ) ioptio = ioptio + 1 431 ! 432 IF( ioptio /= 1 ) & 433 CALL ctl_stop( ' Choose ONE type of light penetration in namelist namtra_qsr', & 434 & ' 2 bands, 3 RGB bands or bio-model light penetration' ) 435 ! 436 IF( ln_qsr_rgb .AND. nn_chldta == 0 ) nqsr = 1 437 IF( ln_qsr_rgb .AND. nn_chldta == 1 ) nqsr = 2 438 IF( ln_qsr_2bd ) nqsr = 3 439 IF( ln_qsr_bio ) nqsr = 4 440 ! 441 IF(lwp) THEN ! Print the choice 442 WRITE(numout,*) 443 IF( nqsr == 1 ) WRITE(numout,*) ' R-G-B light penetration - Constant Chlorophyll' 444 IF( nqsr == 2 ) WRITE(numout,*) ' R-G-B light penetration - Chl data ' 445 IF( nqsr == 3 ) WRITE(numout,*) ' 2 bands light penetration' 446 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' 447 391 ENDIF 448 392 ! 449 ENDIF 450 ! ! ===================================== ! 451 IF( ln_traqsr ) THEN ! Initialisation of Light Penetration ! 452 ! ! ===================================== ! 453 ! 454 xsi0r = 1.e0 / rn_si0 455 xsi1r = 1.e0 / rn_si1 456 ! ! ---------------------------------- ! 457 IF( ln_qsr_rgb ) THEN ! Red-Green-Blue light penetration ! 458 ! ! ---------------------------------- ! 459 ! 460 CALL trc_oce_rgb( rkrgb ) !* tabulated attenuation coef. 461 ! 462 ! !* level of light extinction 463 IF( ln_sco ) THEN ; nksr = jpkm1 464 ELSE ; nksr = trc_oce_ext_lev( r_si2, 0.33e2 ) 465 ENDIF 466 467 IF(lwp) WRITE(numout,*) ' level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m' 468 ! 469 IF( nn_chldta == 1 ) THEN !* Chl data : set sf_chl structure 470 IF(lwp) WRITE(numout,*) 471 IF(lwp) WRITE(numout,*) ' Chlorophyll read in a file' 472 ALLOCATE( sf_chl(1), STAT=ierror ) 473 IF( ierror > 0 ) THEN 474 CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_chl structure' ) ; RETURN 475 ENDIF 476 ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1) ) 477 IF( sn_chl%ln_tint )ALLOCATE( sf_chl(1)%fdta(jpi,jpj,1,2) ) 478 ! ! fill sf_chl with sn_chl and control print 479 CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init', & 480 & 'Solar penetration function of read chlorophyll', 'namtra_qsr' ) 481 ! 482 ELSE !* constant Chl : compute once for all the distribution of light (etot3) 483 IF(lwp) WRITE(numout,*) 484 IF(lwp) WRITE(numout,*) ' Constant Chlorophyll concentration = 0.05' 485 IF( lk_vvl ) THEN ! variable volume 486 IF(lwp) WRITE(numout,*) ' key_vvl: light distribution will be computed at each time step' 487 ELSE ! constant volume: computes one for all 488 IF(lwp) WRITE(numout,*) ' fixed volume: light distribution computed one for all' 489 ! 490 zchl = 0.05 ! constant chlorophyll 491 irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) 492 zekb(:,:) = rkrgb(1,irgb) ! Separation in R-G-B depending of the chlorophyll 493 zekg(:,:) = rkrgb(2,irgb) 494 zekr(:,:) = rkrgb(3,irgb) 495 ! 496 zcoef = ( 1. - rn_abs ) / 3.e0 ! equi-partition in R-G-B 497 ze0(:,:,1) = rn_abs 498 ze1(:,:,1) = zcoef 499 ze2(:,:,1) = zcoef 500 ze3(:,:,1) = zcoef 501 zea(:,:,1) = tmask(:,:,1) ! = ( ze0+ze1+z2+ze3 ) * tmask 502 503 DO jk = 2, nksr+1 504 !CDIR NOVERRCHK 505 DO jj = 1, jpj 506 !CDIR NOVERRCHK 507 DO ji = 1, jpi 508 zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * xsi0r ) 509 zc1 = ze1(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekb(ji,jj) ) 510 zc2 = ze2(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekg(ji,jj) ) 511 zc3 = ze3(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekr(ji,jj) ) 512 ze0(ji,jj,jk) = zc0 513 ze1(ji,jj,jk) = zc1 514 ze2(ji,jj,jk) = zc2 515 ze3(ji,jj,jk) = zc3 516 zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,jk) 517 END DO 518 END DO 519 END DO 520 ! 521 DO jk = 1, nksr 522 etot3(:,:,jk) = r1_rau0_rcp * ( zea(:,:,jk) - zea(:,:,jk+1) ) 523 END DO 524 etot3(:,:,nksr+1:jpk) = 0.e0 ! below 400m set to zero 525 ENDIF 526 ENDIF 527 ! 528 ENDIF 529 ! ! ---------------------------------- ! 530 IF( ln_qsr_2bd ) THEN ! 2 bands light penetration ! 531 ! ! ---------------------------------- ! 532 ! 533 ! ! level of light extinction 534 nksr = trc_oce_ext_lev( rn_si1, 1.e2 ) 535 IF(lwp) THEN 536 WRITE(numout,*) 537 IF(lwp) WRITE(numout,*) ' level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m' 538 ENDIF 539 ! 540 IF( lk_vvl ) THEN ! variable volume 541 IF(lwp) WRITE(numout,*) ' key_vvl: light distribution will be computed at each time step' 542 ELSE ! constant volume: computes one for all 543 zz0 = rn_abs * r1_rau0_rcp 544 zz1 = ( 1. - rn_abs ) * r1_rau0_rcp 545 DO jk = 1, nksr !* solar heat absorbed at T-point computed once for all 546 DO jj = 1, jpj ! top 400 meters 547 DO ji = 1, jpi 548 zc0 = zz0 * EXP( -fsdepw(ji,jj,jk )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk )*xsi1r ) 549 zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r ) 550 etot3(ji,jj,jk) = ( zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1) ) 551 END DO 552 END DO 553 END DO 554 etot3(:,:,nksr+1:jpk) = 0.e0 ! below 400m set to zero 555 ! 556 ENDIF 557 ENDIF 558 ! ! ===================================== ! 559 ELSE ! No light penetration ! 560 ! ! ===================================== ! 561 IF(lwp) THEN 562 WRITE(numout,*) 563 WRITE(numout,*) 'tra_qsr_init : NO solar flux penetration' 564 WRITE(numout,*) '~~~~~~~~~~~~' 565 ENDIF 566 ENDIF 567 ! 568 CALL wrk_dealloc( jpi, jpj, zekb, zekg, zekr ) 569 CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 570 ! 571 IF( nn_timing == 1 ) CALL timing_stop('tra_qsr_init') 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) 410 IF( iom_varid( numror, 'fraqsr_1lev', ldstop = .FALSE. ) > 0 ) THEN 411 CALL iom_get( numror, jpdom_autoglo, 'fraqsr_1lev' , fraqsr_1lev ) 412 ELSE 413 fraqsr_1lev(:,:) = 1._wp ! default : no penetration 414 ENDIF 415 ! 416 IF( nn_timing == 1 ) CALL timing_stop('tra_qsr_init') 572 417 ! 573 418 END SUBROUTINE tra_qsr_init
Note: See TracChangeset
for help on using the changeset viewer.