- Timestamp:
- 2020-05-14T21:46:00+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/TRA/traqsr.F90 (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Property svn:externals
-
old new 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@HEAD sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/TRA/traqsr.F90
r12178 r12928 67 67 68 68 !! * Substitutions 69 # include " vectopt_loop_substitute.h90"69 # include "do_loop_substitute.h90" 70 70 !!---------------------------------------------------------------------- 71 71 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 75 75 CONTAINS 76 76 77 SUBROUTINE tra_qsr( kt )77 SUBROUTINE tra_qsr( kt, Kmm, pts, Krhs ) 78 78 !!---------------------------------------------------------------------- 79 79 !! *** ROUTINE tra_qsr *** … … 87 87 !! I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) ) 88 88 !! The temperature trend associated with the solar radiation penetration 89 !! is given by : zta = 1/e3t dk[ I ] / (r au0*Cp)89 !! is given by : zta = 1/e3t dk[ I ] / (rho0*Cp) 90 90 !! At the bottom, boudary condition for the radiation is no flux : 91 91 !! all heat which has not been absorbed in the above levels is put … … 101 101 !! Morel, A. et Berthon, JF, 1989, Limnol Oceanogr 34(8), 1545-1562 102 102 !!---------------------------------------------------------------------- 103 INTEGER, INTENT(in) :: kt ! ocean time-step 103 INTEGER, INTENT(in ) :: kt ! ocean time-step 104 INTEGER, INTENT(in ) :: Kmm, Krhs ! time level indices 105 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation 104 106 ! 105 107 INTEGER :: ji, jj, jk ! dummy loop indices … … 126 128 IF( l_trdtra ) THEN ! trends diagnostic: save the input temperature trend 127 129 ALLOCATE( ztrdt(jpi,jpj,jpk) ) 128 ztrdt(:,:,:) = tsa(:,:,:,jp_tem)130 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) 129 131 ENDIF 130 132 ! … … 133 135 ! !-----------------------------------! 134 136 IF( kt == nit000 ) THEN !== 1st time step ==! 135 !!gm case neuler not taken into account.... 136 IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN ! read in restart 137 IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 .AND. .NOT.l_1st_euler ) THEN ! read in restart 137 138 IF(lwp) WRITE(numout,*) ' nit000-1 qsr tracer content forcing field read in the restart file' 138 139 z1_2 = 0.5_wp … … 154 155 ! 155 156 DO jk = 1, nksr 156 qsr_hc(:,:,jk) = r1_r au0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) )157 qsr_hc(:,:,jk) = r1_rho0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) ) 157 158 END DO 158 159 ! … … 167 168 DO jk = 1, nksr + 1 168 169 DO jj = 2, jpjm1 ! Separation in R-G-B depending of the surface Chl 169 DO ji = fs_2, fs_jpim1170 DO ji = 2, jpim1 170 171 zchl = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) 171 172 zCtot = 40.6 * zchl**0.459 172 173 zze = 568.2 * zCtot**(-0.746) 173 174 IF( zze > 102. ) zze = 200.0 * zCtot**(-0.293) 174 zpsi = gdepw _n(ji,jj,jk) / zze175 zpsi = gdepw(ji,jj,jk,Kmm) / zze 175 176 ! 176 177 zlogc = LOG( zchl ) … … 195 196 ! 196 197 zcoef = ( 1. - rn_abs ) / 3._wp !* surface equi-partition in R-G-B 197 DO jj = 2, jpjm1 198 DO ji = fs_2, fs_jpim1 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) 204 END DO 198 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) 204 END_2D 205 ! 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 205 226 END DO 206 227 ! 207 DO jk = 2, nksr+1 !* interior equi-partition in R-G-B depending of vertical profile of Chl 208 DO jj = 2, jpjm1 209 DO ji = fs_2, fs_jpim1 210 zchl = MIN( 10. , MAX( 0.03, zchl3d(ji,jj,jk) ) ) 211 irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) 212 zekb(ji,jj) = rkrgb(1,irgb) 213 zekg(ji,jj) = rkrgb(2,irgb) 214 zekr(ji,jj) = rkrgb(3,irgb) 215 END DO 216 END DO 217 218 DO jj = 2, jpjm1 219 DO ji = fs_2, fs_jpim1 220 zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * xsi0r ) 221 zc1 = ze1(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekb(ji,jj) ) 222 zc2 = ze2(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekg(ji,jj) ) 223 zc3 = ze3(ji,jj,jk-1) * EXP( - e3t_n(ji,jj,jk-1) * zekr(ji,jj) ) 224 ze0(ji,jj,jk) = zc0 225 ze1(ji,jj,jk) = zc1 226 ze2(ji,jj,jk) = zc2 227 ze3(ji,jj,jk) = zc3 228 zea(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * wmask(ji,jj,jk) 229 END DO 230 END DO 231 END DO 232 ! 233 DO jk = 1, nksr !* now qsr induced heat content 234 DO jj = 2, jpjm1 235 DO ji = fs_2, fs_jpim1 236 qsr_hc(ji,jj,jk) = r1_rau0_rcp * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) ) 237 END DO 238 END DO 239 END DO 228 DO_3D_00_00( 1, nksr ) 229 qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) ) 230 END_3D 240 231 ! 241 232 DEALLOCATE( zekb , zekg , zekr , ze0 , ze1 , ze2 , ze3 , zea , zchl3d ) … … 243 234 CASE( np_2BD ) !== 2-bands fluxes ==! 244 235 ! 245 zz0 = rn_abs * r1_rau0_rcp ! surface equi-partition in 2-bands 246 zz1 = ( 1. - rn_abs ) * r1_rau0_rcp 247 DO jk = 1, nksr ! solar heat absorbed at T-point in the top 400m 248 DO jj = 2, jpjm1 249 DO ji = fs_2, fs_jpim1 250 zc0 = zz0 * EXP( -gdepw_n(ji,jj,jk )*xsi0r ) + zz1 * EXP( -gdepw_n(ji,jj,jk )*xsi1r ) 251 zc1 = zz0 * EXP( -gdepw_n(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -gdepw_n(ji,jj,jk+1)*xsi1r ) 252 qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) ) 253 END DO 254 END DO 255 END DO 236 zz0 = rn_abs * r1_rho0_rcp ! surface equi-partition in 2-bands 237 zz1 = ( 1. - rn_abs ) * r1_rho0_rcp 238 DO_3D_00_00( 1, nksr ) 239 zc0 = zz0 * EXP( -gdepw(ji,jj,jk ,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk ,Kmm)*xsi1r ) 240 zc1 = zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi1r ) 241 qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0 * wmask(ji,jj,jk) - zc1 * wmask(ji,jj,jk+1) ) 242 END_3D 256 243 ! 257 244 END SELECT 258 245 ! 259 246 ! !-----------------------------! 260 DO jk = 1, nksr ! update to the temp. trend ! 261 DO jj = 2, jpjm1 !-----------------------------! 262 DO ji = fs_2, fs_jpim1 ! vector opt. 263 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) & 264 & + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t_n(ji,jj,jk) 265 END DO 266 END DO 267 END DO 247 DO_3D_00_00( 1, nksr ) 248 pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) & 249 & + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t(ji,jj,jk,Kmm) 250 END_3D 268 251 ! 269 252 ! sea-ice: store the 1st ocean level attenuation coefficient 270 DO jj = 2, jpjm1 271 DO ji = fs_2, fs_jpim1 ! vector opt. 272 IF( qsr(ji,jj) /= 0._wp ) THEN ; fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) ) 273 ELSE ; fraqsr_1lev(ji,jj) = 1._wp 274 ENDIF 275 END DO 276 END DO 253 DO_2D_00_00 254 IF( qsr(ji,jj) /= 0._wp ) THEN ; fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) ) 255 ELSE ; fraqsr_1lev(ji,jj) = 1._wp 256 ENDIF 257 END_2D 277 258 CALL lbc_lnk( 'traqsr', fraqsr_1lev(:,:), 'T', 1._wp ) 278 259 ! … … 281 262 zetot(:,:,nksr+1:jpk) = 0._wp ! below ~400m set to zero 282 263 DO jk = nksr, 1, -1 283 zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * r au0_rcp264 zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * rho0_rcp 284 265 END DO 285 266 CALL iom_put( 'qsr3d', zetot ) ! 3D distribution of shortwave Radiation … … 295 276 ! 296 277 IF( l_trdtra ) THEN ! qsr tracers trends saved for diagnostics 297 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)298 CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt )278 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) 279 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt ) 299 280 DEALLOCATE( ztrdt ) 300 281 ENDIF 301 282 ! ! print mean trends (used for debugging) 302 IF( ln_ctl) CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' qsr - Ta: ', mask1=tmask, clinfo3='tra-ta' )283 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' qsr - Ta: ', mask1=tmask, clinfo3='tra-ta' ) 303 284 ! 304 285 IF( ln_timing ) CALL timing_stop('tra_qsr') … … 336 317 !!---------------------------------------------------------------------- 337 318 ! 338 REWIND( numnam_ref ) ! Namelist namtra_qsr in reference namelist339 319 READ ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901) 340 320 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_qsr in reference namelist' ) 341 321 ! 342 REWIND( numnam_cfg ) ! Namelist namtra_qsr in configuration namelist343 322 READ ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902 ) 344 323 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist' )
Note: See TracChangeset
for help on using the changeset viewer.
