Changeset 13205
- Timestamp:
- 2020-07-02T10:43:08+02:00 (5 years ago)
- Location:
- NEMO/trunk
- Files:
-
- 1 edited
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/TRA/traqsr.F90
r12489 r13205 110 110 REAL(wp) :: zc0 , zc1 , zc2 , zc3 ! - - 111 111 REAL(wp) :: zzc0, zzc1, zzc2, zzc3 ! - - 112 REAL(wp) :: zz0 , zz1 ! - - 113 REAL(wp) :: zCb, zCmax, zze, zpsi, zpsimax, zdelpsi, zCtot, zCze 114 REAL(wp) :: zlogc, zlogc2, zlogc3 115 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zekb, zekg, zekr 116 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt 117 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zetot, zchl3d 112 REAL(wp) :: zz0 , zz1 , ze3t, zlui ! - - 113 REAL(wp) :: zCb, zCmax, zpsi, zpsimax, zrdpsi, zCze 114 REAL(wp) :: zlogc, zlogze, zlogCtot, zlogCze 115 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ze0, ze1, ze2, ze3 116 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt, zetot, ztmp3d 118 117 !!---------------------------------------------------------------------- 119 118 ! … … 160 159 CASE( np_RGB , np_RGBc ) !== R-G-B fluxes ==! 161 160 ! 162 ALLOCATE( ze kb(jpi,jpj) , zekg(jpi,jpj) , zekr (jpi,jpj) ,&163 & ze 0 (jpi,jpj,jpk) , ze1 (jpi,jpj,jpk) , ze2 (jpi,jpj,jpk) ,&164 & z e3 (jpi,jpj,jpk) , zea (jpi,jpj,jpk) , zchl3d(jpi,jpj,jpk) )161 ALLOCATE( ze0 (jpi,jpj) , ze1 (jpi,jpj) , & 162 & ze2 (jpi,jpj) , ze3 (jpi,jpj) , & 163 & ztmp3d(jpi,jpj,nksr + 1) ) 165 164 ! 166 165 IF( nqsr == np_RGBc ) THEN !* Variable Chlorophyll 167 166 CALL fld_read( kt, 1, sf_chl ) ! Read Chl data and provides it at the current time step 167 ! 168 ! Separation in R-G-B depending on the surface Chl 169 ! perform and store as many of the 2D calculations as possible 170 ! before the 3D loop (use the temporary 2D arrays to replace the 171 ! most expensive calculations) 172 ! 173 DO_2D_00_00 174 ! zlogc = log(zchl) 175 zlogc = LOG ( MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) ) 176 ! zc1 : log(zCze) = log (1.12 * zchl**0.803) 177 zc1 = 0.113328685307 + 0.803 * zlogc 178 ! zc2 : log(zCtot) = log(40.6 * zchl**0.459) 179 zc2 = 3.703768066608 + 0.459 * zlogc 180 ! zc3 : log(zze) = log(568.2 * zCtot**(-0.746)) 181 zc3 = 6.34247346942 - 0.746 * zc2 182 ! IF( log(zze) > log(102.) ) log(zze) = log(200.0 * zCtot**(-0.293)) 183 IF( zc3 > 4.62497281328 ) zc3 = 5.298317366548 - 0.293 * zc2 184 ! 185 ze0(ji,jj) = zlogc ! ze0 = log(zchl) 186 ze1(ji,jj) = EXP( zc1 ) ! ze1 = zCze 187 ze2(ji,jj) = 1._wp / ( 0.710 + zlogc * ( 0.159 + zlogc * 0.021 ) ) ! ze2 = 1/zdelpsi 188 ze3(ji,jj) = EXP( - zc3 ) ! ze3 = 1/zze 189 END_2D 190 191 ! 192 DO_3D_00_00 ( 1, nksr + 1 ) 193 ! zchl = ALOG( ze0(ji,jj) ) 194 zlogc = ze0(ji,jj) 195 ! 196 zCb = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) ) 197 zCmax = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 ) 198 zpsimax = 0.6 - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) ) 199 ! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 ) 200 ! 201 zCze = ze1(ji,jj) 202 zrdpsi = ze2(ji,jj) ! 1/zdelpsi 203 zpsi = ze3(ji,jj) * gdepw(ji,jj,jk,Kmm) ! gdepw/zze 204 ! 205 ! NB. make sure zchl value is such that: zchl = MIN( 10. , MAX( 0.03, zchl ) ) 206 zchl = MIN( 10. , MAX( 0.03, zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) ) ) 207 ! Convert chlorophyll value to attenuation coefficient look-up table index 208 ztmp3d(ji,jj,jk) = 41 + 20.*LOG10(zchl) + 1.e-15 209 END_3D 210 ELSE !* constant chlorophyll 211 zchl = 0.05 212 ! NB. make sure constant value is such that: 213 zchl = MIN( 10. , MAX( 0.03, zchl ) ) 214 ! Convert chlorophyll value to attenuation coefficient look-up table index 215 zlui = 41 + 20.*LOG10(zchl) + 1.e-15 168 216 DO jk = 1, nksr + 1 169 DO jj = 2, jpjm1 ! Separation in R-G-B depending of the surface Chl 170 DO ji = 2, jpim1 171 zchl = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) 172 zCtot = 40.6 * zchl**0.459 173 zze = 568.2 * zCtot**(-0.746) 174 IF( zze > 102. ) zze = 200.0 * zCtot**(-0.293) 175 zpsi = gdepw(ji,jj,jk,Kmm) / zze 176 ! 177 zlogc = LOG( zchl ) 178 zlogc2 = zlogc * zlogc 179 zlogc3 = zlogc * zlogc * zlogc 180 zCb = 0.768 + 0.087 * zlogc - 0.179 * zlogc2 - 0.025 * zlogc3 181 zCmax = 0.299 - 0.289 * zlogc + 0.579 * zlogc2 182 zpsimax = 0.6 - 0.640 * zlogc + 0.021 * zlogc2 + 0.115 * zlogc3 183 zdelpsi = 0.710 + 0.159 * zlogc + 0.021 * zlogc2 184 zCze = 1.12 * (zchl)**0.803 185 ! 186 zchl3d(ji,jj,jk) = zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) / zdelpsi )**2 ) ) 187 END DO 188 ! 189 END DO 217 ztmp3d(:,:,jk) = zlui 190 218 END DO 191 ELSE !* constant chrlorophyll192 DO jk = 1, nksr + 1193 zchl3d(:,:,jk) = 0.05194 ENDDO195 219 ENDIF 196 220 ! 197 221 zcoef = ( 1. - rn_abs ) / 3._wp !* surface equi-partition in R-G-B 198 222 DO_2D_00_00 199 ze0(ji,jj,1) = rn_abs * qsr(ji,jj) 200 ze1(ji,jj,1) = zcoef * qsr(ji,jj) 201 ze2(ji,jj,1) = zcoef * qsr(ji,jj) 202 ze3(ji,jj,1) = zcoef * qsr(ji,jj) 203 zea(ji,jj,1) = qsr(ji,jj) 223 ze0(ji,jj) = rn_abs * qsr(ji,jj) 224 ze1(ji,jj) = zcoef * qsr(ji,jj) 225 ze2(ji,jj) = zcoef * qsr(ji,jj) 226 ze3(ji,jj) = zcoef * qsr(ji,jj) 227 ! store the surface SW radiation; re-use the surface ztmp3d array 228 ! since the surface attenuation coefficient is not used 229 ztmp3d(ji,jj,1) = qsr(ji,jj) 204 230 END_2D 205 231 ! 206 DO jk = 2, nksr+1 !* interior equi-partition in R-G-B depending of vertical profile of Chl 207 DO_2D_00_00 208 zchl = MIN( 10. , MAX( 0.03, zchl3d(ji,jj,jk) ) ) 209 irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) 210 zekb(ji,jj) = rkrgb(1,irgb) 211 zekg(ji,jj) = rkrgb(2,irgb) 212 zekr(ji,jj) = rkrgb(3,irgb) 213 END_2D 214 215 DO_2D_00_00 216 zc0 = ze0(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * xsi0r ) 217 zc1 = ze1(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * zekb(ji,jj) ) 218 zc2 = ze2(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * zekg(ji,jj) ) 219 zc3 = ze3(ji,jj,jk-1) * EXP( - e3t(ji,jj,jk-1,Kmm) * zekr(ji,jj) ) 220 ze0(ji,jj,jk) = zc0 221 ze1(ji,jj,jk) = zc1 222 ze2(ji,jj,jk) = zc2 223 ze3(ji,jj,jk) = zc3 224 zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk) 225 END_2D 226 END DO 232 !* interior equi-partition in R-G-B depending on vertical profile of Chl 233 DO_3D_00_00 ( 2, nksr + 1 ) 234 ze3t = e3t(ji,jj,jk-1,Kmm) 235 irgb = NINT( ztmp3d(ji,jj,jk) ) 236 zc0 = ze0(ji,jj) * EXP( - ze3t * xsi0r ) 237 zc1 = ze1(ji,jj) * EXP( - ze3t * rkrgb(1,irgb) ) 238 zc2 = ze2(ji,jj) * EXP( - ze3t * rkrgb(2,irgb) ) 239 zc3 = ze3(ji,jj) * EXP( - ze3t * rkrgb(3,irgb) ) 240 ze0(ji,jj) = zc0 241 ze1(ji,jj) = zc1 242 ze2(ji,jj) = zc2 243 ze3(ji,jj) = zc3 244 ztmp3d(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk) 245 END_3D 227 246 ! 228 247 DO_3D_00_00( 1, nksr ) 229 qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( z ea(ji,jj,jk) - zea(ji,jj,jk+1) )248 qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( ztmp3d(ji,jj,jk) - ztmp3d(ji,jj,jk+1) ) 230 249 END_3D 231 250 ! 232 DEALLOCATE( ze kb , zekg , zekr , ze0 , ze1 , ze2 , ze3 , zea , zchl3d )251 DEALLOCATE( ze0 , ze1 , ze2 , ze3 , ztmp3d ) 233 252 ! 234 253 CASE( np_2BD ) !== 2-bands fluxes ==!
Note: See TracChangeset
for help on using the changeset viewer.
