- Timestamp:
- 2020-06-10T13:13:39+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
r11442 r13088 9 9 !! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module 10 10 !! - ! 2005-11 (G. Madec) zco, zps, sco coordinate 11 !! 3.2 ! 2009-04 (G. Madec & NEMO team) 12 !! 3.4 ! 2012-05 (C. Rousset) store attenuation coef for use in ice model 11 !! 3.2 ! 2009-04 (G. Madec & NEMO team) 12 !! 3.4 ! 2012-05 (C. Rousset) store attenuation coef for use in ice model 13 13 !! 3.6 ! 2015-12 (O. Aumont, J. Jouanno, C. Ethe) use vertical profile of chlorophyll 14 14 !!---------------------------------------------------------------------- … … 43 43 ! !!* Namelist namtra_qsr: penetrative solar radiation 44 44 LOGICAL , PUBLIC :: ln_traqsr !: light absorption (qsr) flag 45 LOGICAL , PUBLIC :: ln_qsr_rgb !: Red-Green-Blue light absorption flag 45 LOGICAL , PUBLIC :: ln_qsr_rgb !: Red-Green-Blue light absorption flag 46 46 LOGICAL , PUBLIC :: ln_qsr_2bd !: 2 band light absorption flag 47 47 LOGICAL , PUBLIC :: ln_qsr_bio !: bio-model light absorption flag … … 51 51 REAL(wp), PUBLIC :: rn_si0 !: very near surface depth of extinction (RGB & 2 bands) 52 52 REAL(wp), PUBLIC :: rn_si1 !: deepest depth of extinction (water type I) (2 bands) 53 53 54 54 ! Module variables 55 55 REAL(wp), ALLOCATABLE :: xsi0r(:,:) !: inverse of rn_si0 … … 80 80 !! Considering the 2 wavebands case: 81 81 !! I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) ) 82 !! The temperature trend associated with the solar radiation penetration 82 !! The temperature trend associated with the solar radiation penetration 83 83 !! is given by : zta = 1/e3t dk[ I ] / (rau0*Cp) 84 84 !! At the bottom, boudary condition for the radiation is no flux : … … 86 86 !! in the last ocean level. 87 87 !! In z-coordinate case, the computation is only done down to the 88 !! level where I(k) < 1.e-15 W/m2. In addition, the coefficients 88 !! level where I(k) < 1.e-15 W/m2. In addition, the coefficients 89 89 !! used for the computation are calculated one for once as they 90 90 !! depends on k only. … … 106 106 REAL(wp) :: zz0, zz1, z1_e3t ! - - 107 107 REAL(wp) :: zCb, zCmax, zze, zpsi, zpsimax, zdelpsi, zCtot, zCze 108 REAL(wp) :: zlogc, zlogc2, zlogc3 108 REAL(wp) :: zlogc, zlogc2, zlogc3 109 109 REAL(wp), POINTER, DIMENSION(:,: ) :: zekb, zekg, zekr 110 110 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt, zchl3d … … 113 113 IF( nn_timing == 1 ) CALL timing_start('tra_qsr') 114 114 ! 115 CALL wrk_alloc( jpi, jpj, zekb, zekg, zekr ) 116 CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d ) 115 CALL wrk_alloc( jpi, jpj, zekb, zekg, zekr ) 116 CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d ) 117 117 ! 118 118 IF( kt == nit000 ) THEN … … 124 124 125 125 IF( l_trdtra ) THEN ! Save ta and sa trends 126 CALL wrk_alloc( jpi, jpj, jpk, ztrdt ) 126 CALL wrk_alloc( jpi, jpj, jpk, ztrdt ) 127 127 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 128 128 ENDIF … … 150 150 ! Compute now qsr tracer content field 151 151 ! ************************************ 152 152 153 153 ! ! ============================================== ! 154 154 IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN ! bio-model fluxes : all vertical coordinates ! … … 159 159 ! Add to the general trend 160 160 DO jk = 1, jpkm1 161 DO jj = 2, jpjm1 161 DO jj = 2, jpjm1 162 162 DO ji = fs_2, fs_jpim1 ! vector opt. 163 163 z1_e3t = zfact / fse3t(ji,jj,jk) … … 180 180 ENDIF 181 181 ! ! ============================================== ! 182 ELSE ! Ocean alone : 182 ELSE ! Ocean alone : 183 183 ! ! ============================================== ! 184 184 ! 185 ! 185 ! 186 #if defined key_traldf_c2d || key_traldf_c3d 186 187 IF( ln_stopack .AND. ( nn_spp_qsi0 > 0 ) ) THEN 187 xsi0r = rn_si0 188 CALL spp_gen(kt, xsi0r, nn_spp_qsi0, rn_qsi0_sd, jk_spp_qsi0 ) 189 xsi0r = 1.e0 / xsi0r 190 ENDIF 188 xsi0r = rn_si0 189 CALL spp_gen(kt, xsi0r, nn_spp_qsi0, rn_qsi0_sd, jk_spp_qsi0 ) 190 xsi0r = 1.e0 / xsi0r 191 ENDIF 192 #else 193 IF ( ln_stopack .AND. nn_spp_qsi0 > 0 ) & 194 & CALL ctl_stop( 'tra_qsr: parameter perturbation will only work with '// & 195 'key_traldf_c2d or key_traldf_c3d') 196 #endif 197 191 198 ! ! ------------------------- ! 192 199 IF( ln_qsr_rgb) THEN ! R-G-B light penetration ! … … 199 206 CALL fld_read( kt, 1, sf_chl ) ! Read Chl data and provides it at the current time step 200 207 DO jk = 1, nksr + 1 201 zchl3d(:,:,jk) = sf_chl(1)%fnow(:,:,1) 208 zchl3d(:,:,jk) = sf_chl(1)%fnow(:,:,1) 202 209 ENDDO 203 210 ! … … 220 227 zpsimax = 0.6 - 0.640 * zlogc + 0.021 * zlogc2 + 0.115 * zlogc3 221 228 zdelpsi = 0.710 + 0.159 * zlogc + 0.021 * zlogc2 222 zCze = 1.12 * (zchl)**0.803 229 zCze = 1.12 * (zchl)**0.803 223 230 DO jk = 1, nksr + 1 224 231 zpsi = fsdept(ji,jj,jk) / zze … … 231 238 ELSE !* Variable ocean volume but constant chrlorophyll 232 239 DO jk = 1, nksr + 1 233 zchl3d(:,:,jk) = 0.05 240 zchl3d(:,:,jk) = 0.05 234 241 ENDDO 235 242 ENDIF … … 256 263 !CDIR NOVERRCHK 257 264 DO jj = 1, jpj 258 !CDIR NOVERRCHK 265 !CDIR NOVERRCHK 259 266 DO ji = 1, jpi 260 267 zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r(ji,jj) ) … … 289 296 END DO 290 297 END DO 291 ! 298 ! 292 299 DO jj = 1, jpj 293 300 DO ji = 1, jpi … … 296 303 zc2 = zcoef * EXP( - fse3t(ji,jj,1) * zekg(ji,jj) ) 297 304 zc3 = zcoef * EXP( - fse3t(ji,jj,1) * zekr(ji,jj) ) 298 fraqsr_1lev(ji,jj) = 1.0 - ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,2) 305 fraqsr_1lev(ji,jj) = 1.0 - ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,2) 299 306 END DO 300 307 END DO … … 321 328 zz0 = rn_abs * r1_rau0_rcp 322 329 zz1 = ( 1. - rn_abs ) * r1_rau0_rcp 323 DO jk = 1, nksr ! solar heat absorbed at T-point in the top 400m 330 DO jk = 1, nksr ! solar heat absorbed at T-point in the top 400m 324 331 DO jj = 1, jpj 325 332 DO ji = 1, jpi 326 333 zc0 = zz0 * EXP( -fsdepw(ji,jj,jk )*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk )*xsi1r ) 327 334 zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r ) 328 qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) ) 335 qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) ) 329 336 END DO 330 337 END DO … … 344 351 DO jj = 2, jpjm1 345 352 DO ji = fs_2, fs_jpim1 ! vector opt. 346 ! (ISF) no light penetration below the ice shelves 353 ! (ISF) no light penetration below the ice shelves 347 354 qsr_hc(ji,jj,jk) = etot3(ji,jj,jk) * qsr(ji,jj) * tmask(ji,jj,1) 348 355 END DO … … 360 367 ! Add to the general trend 361 368 DO jk = 1, nksr 362 DO jj = 2, jpjm1 369 DO jj = 2, jpjm1 363 370 DO ji = fs_2, fs_jpim1 ! vector opt. 364 371 z1_e3t = zfact / fse3t(ji,jj,jk) … … 376 383 & 'at it= ', kt,' date= ', ndastp 377 384 IF(lwp) WRITE(numout,*) '~~~~' 378 IF(nn_timing == 2) CALL timing_start('iom_rstput') 385 IF(nn_timing == 2) CALL timing_start('iom_rstput') 379 386 CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b' , qsr_hc ) 380 CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev ) ! default definition in sbcssm 387 CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev ) ! default definition in sbcssm 381 388 IF(nn_timing == 2) CALL timing_stop('iom_rstput') 382 389 ! … … 386 393 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 387 394 CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt ) 388 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt ) 395 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt ) 389 396 ENDIF 390 397 ! ! print mean trends (used for debugging) 391 398 IF(ln_ctl) CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' qsr - Ta: ', mask1=tmask, clinfo3='tra-ta' ) 392 399 ! 393 CALL wrk_dealloc( jpi, jpj, zekb, zekg, zekr ) 394 CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d ) 400 CALL wrk_dealloc( jpi, jpj, zekb, zekg, zekr ) 401 CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d ) 395 402 ! 396 403 IF( nn_timing == 1 ) CALL timing_stop('tra_qsr') … … 408 415 !! from two length scale of penetration (rn_si0,rn_si1) and a ratio 409 416 !! (rn_abs). These parameters are read in the namtra_qsr namelist. The 410 !! default values correspond to clear water (type I in Jerlov' 417 !! default values correspond to clear water (type I in Jerlov' 411 418 !! (1968) classification. 412 419 !! called by tra_qsr at the first timestep (nit000) … … 435 442 IF( nn_timing == 1 ) CALL timing_start('tra_qsr_init') 436 443 ! 437 CALL wrk_alloc( jpi, jpj, zekb, zekg, zekr ) 438 CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 444 CALL wrk_alloc( jpi, jpj, zekb, zekg, zekr ) 445 CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 439 446 ! 440 447 … … 465 472 466 473 IF( ln_traqsr ) THEN ! control consistency 467 ! 474 ! 468 475 IF( .NOT.lk_qsr_bio .AND. ln_qsr_bio ) THEN 469 476 CALL ctl_warn( 'No bio model : force ln_qsr_bio = FALSE ' ) … … 480 487 & ' 2 bands, 3 RGB bands or bio-model light penetration' ) 481 488 ! 482 IF( ln_qsr_rgb .AND. nn_chldta == 0 ) nqsr = 1 489 IF( ln_qsr_rgb .AND. nn_chldta == 0 ) nqsr = 1 483 490 IF( ln_qsr_rgb .AND. nn_chldta == 1 ) nqsr = 2 484 491 IF( ln_qsr_rgb .AND. nn_chldta == 2 ) nqsr = 3 … … 497 504 ENDIF 498 505 ! ! ===================================== ! 499 IF( ln_traqsr ) THEN ! Initialisation of Light Penetration ! 506 IF( ln_traqsr ) THEN ! Initialisation of Light Penetration ! 500 507 ! ! ===================================== ! 501 508 ! … … 539 546 zchl = 0.05 ! constant chlorophyll 540 547 irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) 541 zekb(:,:) = rkrgb(1,irgb) ! Separation in R-G-B depending of the chlorophyll 548 zekb(:,:) = rkrgb(1,irgb) ! Separation in R-G-B depending of the chlorophyll 542 549 zekg(:,:) = rkrgb(2,irgb) 543 550 zekr(:,:) = rkrgb(3,irgb) … … 546 553 ze0(:,:,1) = rn_abs 547 554 ze1(:,:,1) = zcoef 548 ze2(:,:,1) = zcoef 555 ze2(:,:,1) = zcoef 549 556 ze3(:,:,1) = zcoef 550 557 zea(:,:,1) = tmask(:,:,1) ! = ( ze0+ze1+z2+ze3 ) * tmask 551 558 552 559 DO jk = 2, nksr+1 553 560 !CDIR NOVERRCHK 554 561 DO jj = 1, jpj 555 !CDIR NOVERRCHK 562 !CDIR NOVERRCHK 556 563 DO ji = 1, jpi 557 564 zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * xsi0r(ji,jj) ) … … 566 573 END DO 567 574 END DO 568 END DO 575 END DO 569 576 ! 570 577 DO jk = 1, nksr … … 598 605 zc0 = zz0 * EXP( -fsdepw(ji,jj,jk )*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk )*xsi1r ) 599 606 zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r ) 600 etot3(ji,jj,jk) = ( zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1) ) * tmask(ji,jj,1) 607 etot3(ji,jj,jk) = ( zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1) ) * tmask(ji,jj,1) 601 608 END DO 602 609 END DO … … 607 614 ENDIF 608 615 ! ! ===================================== ! 609 ELSE ! No light penetration ! 616 ELSE ! No light penetration ! 610 617 ! ! ===================================== ! 611 618 IF(lwp) THEN … … 625 632 ENDIF 626 633 ! 627 CALL wrk_dealloc( jpi, jpj, zekb, zekg, zekr ) 628 CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 634 CALL wrk_dealloc( jpi, jpj, zekb, zekg, zekr ) 635 CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 629 636 ! 630 637 IF( nn_timing == 1 ) CALL timing_stop('tra_qsr_init')
Note: See TracChangeset
for help on using the changeset viewer.
