Changeset 12956
- Timestamp:
- 2020-05-20T16:24:20+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12953_ENHANCE-10_acc_fix_traqsr/src/OCE/TRA/traqsr.F90
r12489 r12956 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 ze0(:,:) = LOG ( MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(:,:,1) ) ) ) ! ze0 = log(zchl) 174 ze1(:,:) = 0.113328685307 + 0.803 * ze0(:,:) ! ze1 : log(zCze) = log (1.12 * zchl**0.803) 175 ze2(:,:) = 3.703768066608 + 0.459 * ze0(:,:) ! ze2 : log(zCtot) = log(40.6 * zchl**0.459) 176 ze3(:,:) = 6.34247346942 - 0.746 * ze2(:,:) ! ze3 : log(zze) = log(568.2 * zCtot**(-0.746)) 177 WHERE( ze3 > 4.62497281328 ) ze3 = 5.298317366548 - 0.293 * ze2 ! IF( log(zze) > log(102.) ) log(zze) = 178 ! ! log(200.0 * zCtot**(-0.293)) 179 ze1(:,:) = EXP(ze1(:,:)) ! ze1 = zCze 180 ze2(:,:) = 1._wp / ( 0.710 + ze0(:,:) * ( 0.159 + ze0(:,:) * 0.021 ) ) ! ze2 = 1/zdelpsi 181 ze3(:,:) = 1._wp / EXP(ze3(:,:)) ! ze3 = 1/zze 182 ! 183 DO_3D_00_00 ( 1, nksr + 1 ) 184 ! zchl = ALOG( ze0(ji,jj) ) 185 zlogc = ze0(ji,jj) 186 ! 187 zCb = 0.768 + zlogc * ( 0.087 - zlogc * ( 0.179 + zlogc * 0.025 ) ) 188 zCmax = 0.299 - zlogc * ( 0.289 - zlogc * 0.579 ) 189 zpsimax = 0.6 - zlogc * ( 0.640 - zlogc * ( 0.021 + zlogc * 0.115 ) ) 190 ! zdelpsi = 0.710 + zlogc * ( 0.159 + zlogc * 0.021 ) 191 ! 192 zCze = ze1(ji,jj) 193 zrdpsi = ze2(ji,jj) ! 1/zdelpsi 194 zpsi = ze3(ji,jj) * gdepw(ji,jj,jk,Kmm) ! gdepw/zze 195 ! 196 ! NB. make sure zchl value is such that: zchl = MIN( 10. , MAX( 0.03, zchl ) ) 197 zchl = MIN( 10. , MAX( 0.03, zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) * zrdpsi )**2 ) ) ) ) 198 ! Convert chlorophyll value to attenuation coefficient look-up table index 199 ztmp3d(ji,jj,jk) = 41 + 20.*LOG10(zchl) + 1.e-15 200 END_3D 201 ELSE !* constant chlorophyll 202 zchl = 0.05 203 ! NB. make sure constant value is such that: 204 zchl = MIN( 10. , MAX( 0.03, zchl ) ) 205 ! Convert chlorophyll value to attenuation coefficient look-up table index 206 zlui = 41 + 20.*LOG10(zchl) + 1.e-15 168 207 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 208 ztmp3d(:,:,jk) = zlui 190 209 END DO 191 ELSE !* constant chrlorophyll192 DO jk = 1, nksr + 1193 zchl3d(:,:,jk) = 0.05194 ENDDO195 210 ENDIF 196 211 ! 197 212 zcoef = ( 1. - rn_abs ) / 3._wp !* surface equi-partition in R-G-B 198 213 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) 214 ze0(ji,jj) = rn_abs * qsr(ji,jj) 215 ze1(ji,jj) = zcoef * qsr(ji,jj) 216 ze2(ji,jj) = zcoef * qsr(ji,jj) 217 ze3(ji,jj) = zcoef * qsr(ji,jj) 218 ! store the surface SW radiation; re-use the surface ztmp3d array 219 ! since the surface attenuation coefficient is not used 220 ztmp3d(ji,jj,1) = qsr(ji,jj) 204 221 END_2D 205 222 ! 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 223 !* interior equi-partition in R-G-B depending on vertical profile of Chl 224 DO_3D_00_00 ( 2, nksr + 1 ) 225 ze3t = e3t(ji,jj,jk-1,Kmm) 226 irgb = NINT( ztmp3d(ji,jj,jk) ) 227 zc0 = ze0(ji,jj) * EXP( - ze3t * xsi0r ) 228 zc1 = ze1(ji,jj) * EXP( - ze3t * rkrgb(1,irgb) ) 229 zc2 = ze2(ji,jj) * EXP( - ze3t * rkrgb(2,irgb) ) 230 zc3 = ze3(ji,jj) * EXP( - ze3t * rkrgb(3,irgb) ) 231 ze0(ji,jj) = zc0 232 ze1(ji,jj) = zc1 233 ze2(ji,jj) = zc2 234 ze3(ji,jj) = zc3 235 ztmp3d(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk) 236 END_3D 227 237 ! 228 238 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) )239 qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( ztmp3d(ji,jj,jk) - ztmp3d(ji,jj,jk+1) ) 230 240 END_3D 231 241 ! 232 DEALLOCATE( ze kb , zekg , zekr , ze0 , ze1 , ze2 , ze3 , zea , zchl3d )242 DEALLOCATE( ze0 , ze1 , ze2 , ze3 , ztmp3d ) 233 243 ! 234 244 CASE( np_2BD ) !== 2-bands fluxes ==!
Note: See TracChangeset
for help on using the changeset viewer.
