Changeset 15603 for branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
- Timestamp:
- 2021-12-16T10:39:55+01:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
r12555 r15603 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 !!---------------------------------------------------------------------- … … 33 33 USE wrk_nemo ! Memory Allocation 34 34 USE timing ! Timing 35 USE stopack 35 36 36 37 IMPLICIT NONE … … 42 43 ! !!* Namelist namtra_qsr: penetrative solar radiation 43 44 LOGICAL , PUBLIC :: ln_traqsr !: light absorption (qsr) flag 44 LOGICAL , PUBLIC :: ln_qsr_rgb !: Red-Green-Blue light absorption flag 45 LOGICAL , PUBLIC :: ln_qsr_rgb !: Red-Green-Blue light absorption flag 45 46 LOGICAL , PUBLIC :: ln_qsr_2bd !: 2 band light absorption flag 46 47 LOGICAL , PUBLIC :: ln_qsr_bio !: bio-model light absorption flag … … 50 51 REAL(wp), PUBLIC :: rn_si0 !: very near surface depth of extinction (RGB & 2 bands) 51 52 REAL(wp), PUBLIC :: rn_si1 !: deepest depth of extinction (water type I) (2 bands) 52 53 53 54 ! Module variables 54 REAL(wp) :: xsi0r!: inverse of rn_si055 REAL(wp), ALLOCATABLE :: xsi0r(:,:) !: inverse of rn_si0 55 56 REAL(wp) :: xsi1r !: inverse of rn_si1 56 57 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_chl ! structure of input Chl (file informations, fields read) … … 79 80 !! Considering the 2 wavebands case: 80 81 !! I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) ) 81 !! The temperature trend associated with the solar radiation penetration 82 !! The temperature trend associated with the solar radiation penetration 82 83 !! is given by : zta = 1/e3t dk[ I ] / (rau0*Cp) 83 84 !! At the bottom, boudary condition for the radiation is no flux : … … 85 86 !! in the last ocean level. 86 87 !! 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 !! level where I(k) < 1.e-15 W/m2. In addition, the coefficients 88 89 !! used for the computation are calculated one for once as they 89 90 !! depends on k only. … … 105 106 REAL(wp) :: zz0, zz1, z1_e3t ! - - 106 107 REAL(wp) :: zCb, zCmax, zze, zpsi, zpsimax, zdelpsi, zCtot, zCze 107 REAL(wp) :: zlogc, zlogc2, zlogc3 108 REAL(wp) :: zlogc, zlogc2, zlogc3 108 109 REAL(wp), POINTER, DIMENSION(:,: ) :: zekb, zekg, zekr 109 110 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt, zchl3d … … 112 113 IF( nn_timing == 1 ) CALL timing_start('tra_qsr') 113 114 ! 114 CALL wrk_alloc( jpi, jpj, zekb, zekg, zekr ) 115 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 ) 116 117 ! 117 118 IF( kt == nit000 ) THEN … … 124 125 125 126 IF( l_trdtra ) THEN ! Save ta and sa trends 126 CALL wrk_alloc( jpi, jpj, jpk, ztrdt ) 127 CALL wrk_alloc( jpi, jpj, jpk, ztrdt ) 127 128 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 128 129 ENDIF … … 153 154 ! Compute now qsr tracer content field 154 155 ! ************************************ 155 156 156 157 ! ! ============================================== ! 157 158 IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN ! bio-model fluxes : all vertical coordinates ! … … 162 163 ! Add to the general trend 163 164 DO jk = 1, jpkm1 164 DO jj = 2, jpjm1 165 DO jj = 2, jpjm1 165 166 DO ji = fs_2, fs_jpim1 ! vector opt. 166 167 z1_e3t = zfact / fse3t(ji,jj,jk) … … 183 184 ENDIF 184 185 ! ! ============================================== ! 185 ELSE ! Ocean alone : 186 ELSE ! Ocean alone : 186 187 ! ! ============================================== ! 187 188 ! 188 ! ! ------------------------- ! 189 ! 190 #if defined key_traldf_c2d || key_traldf_c3d 191 IF( ln_stopack .AND. ( nn_spp_qsi0 > 0 ) ) THEN 192 xsi0r = rn_si0 193 CALL spp_gen(kt, xsi0r, nn_spp_qsi0, rn_qsi0_sd, jk_spp_qsi0 ) 194 xsi0r = 1.e0 / xsi0r 195 ENDIF 196 #else 197 IF ( ln_stopack .AND. nn_spp_qsi0 > 0 ) & 198 & CALL ctl_stop( 'tra_qsr: parameter perturbation will only work with '// & 199 'key_traldf_c2d or key_traldf_c3d') 200 #endif 201 202 ! ! ------------------------- ! 189 203 IF( ln_qsr_rgb) THEN ! R-G-B light penetration ! 190 204 ! ! ------------------------- ! … … 196 210 CALL fld_read( kt, 1, sf_chl ) ! Read Chl data and provides it at the current time step 197 211 DO jk = 1, nksr + 1 198 zchl3d(:,:,jk) = sf_chl(1)%fnow(:,:,1) 212 zchl3d(:,:,jk) = sf_chl(1)%fnow(:,:,1) 199 213 ENDDO 200 214 ! … … 217 231 zpsimax = 0.6 - 0.640 * zlogc + 0.021 * zlogc2 + 0.115 * zlogc3 218 232 zdelpsi = 0.710 + 0.159 * zlogc + 0.021 * zlogc2 219 zCze = 1.12 * (zchl)**0.803 233 zCze = 1.12 * (zchl)**0.803 220 234 DO jk = 1, nksr + 1 221 235 zpsi = fsdept(ji,jj,jk) / zze … … 228 242 ELSE !* Variable ocean volume but constant chrlorophyll 229 243 DO jk = 1, nksr + 1 230 zchl3d(:,:,jk) = 0.05 244 zchl3d(:,:,jk) = 0.05 231 245 ENDDO 232 246 ENDIF … … 253 267 !CDIR NOVERRCHK 254 268 DO jj = 1, jpj 255 !CDIR NOVERRCHK 256 DO ji = 1, jpi 257 zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r 269 !CDIR NOVERRCHK 270 DO ji = 1, jpi 271 zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r(ji,jj) ) 258 272 zc1 = ze1(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekb(ji,jj) ) 259 273 zc2 = ze2(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * zekg(ji,jj) ) … … 286 300 END DO 287 301 END DO 288 ! 302 ! 289 303 DO jj = 1, jpj 290 304 DO ji = 1, jpi 291 zc0 = rn_abs * EXP( - fse3t(ji,jj,1) * xsi0r 305 zc0 = rn_abs * EXP( - fse3t(ji,jj,1) * xsi0r(ji,jj) ) 292 306 zc1 = zcoef * EXP( - fse3t(ji,jj,1) * zekb(ji,jj) ) 293 307 zc2 = zcoef * EXP( - fse3t(ji,jj,1) * zekg(ji,jj) ) 294 308 zc3 = zcoef * EXP( - fse3t(ji,jj,1) * zekr(ji,jj) ) 295 fraqsr_1lev(ji,jj) = 1.0 - ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,2) 309 fraqsr_1lev(ji,jj) = 1.0 - ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,2) 296 310 END DO 297 311 END DO … … 314 328 ! ! ------------------------- ! 315 329 ! 316 IF( lk_vvl ) THEN !* variable volume 330 IF( lk_vvl .OR. ( ln_stopack .AND. ( nn_spp_qsi0 > 0 ) ) ) THEN !* variable volume 331 317 332 zz0 = rn_abs * r1_rau0_rcp 318 333 zz1 = ( 1. - rn_abs ) * r1_rau0_rcp 319 DO jk = 1, nksr ! solar heat absorbed at T-point in the top 400m 334 DO jk = 1, nksr ! solar heat absorbed at T-point in the top 400m 320 335 DO jj = 1, jpj 321 336 DO ji = 1, jpi 322 zc0 = zz0 * EXP( -fsdepw(ji,jj,jk )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk )*xsi1r )323 zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r )324 qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) ) 337 zc0 = zz0 * EXP( -fsdepw(ji,jj,jk )*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk )*xsi1r ) 338 zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r ) 339 qsr_hc(ji,jj,jk) = qsr(ji,jj) * ( zc0*tmask(ji,jj,jk) - zc1*tmask(ji,jj,jk+1) ) 325 340 END DO 326 341 END DO … … 330 345 DO jj = 1, jpj 331 346 DO ji = 1, jpi 332 zc0 = zz0 * EXP( -fsdepw(ji,jj,1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,1)*xsi1r )333 zc1 = zz0 * EXP( -fsdepw(ji,jj,2)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,2)*xsi1r )347 zc0 = zz0 * EXP( -fsdepw(ji,jj,1)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,1)*xsi1r ) 348 zc1 = zz0 * EXP( -fsdepw(ji,jj,2)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,2)*xsi1r ) 334 349 fraqsr_1lev(ji,jj) = ( zc0*tmask(ji,jj,1) - zc1*tmask(ji,jj,2) ) / r1_rau0_rcp 335 350 END DO … … 340 355 DO jj = 2, jpjm1 341 356 DO ji = fs_2, fs_jpim1 ! vector opt. 342 ! (ISF) no light penetration below the ice shelves 357 ! (ISF) no light penetration below the ice shelves 343 358 qsr_hc(ji,jj,jk) = etot3(ji,jj,jk) * qsr(ji,jj) * tmask(ji,jj,1) 344 359 END DO … … 356 371 ! Add to the general trend 357 372 DO jk = 1, nksr 358 DO jj = 2, jpjm1 373 DO jj = 2, jpjm1 359 374 DO ji = fs_2, fs_jpim1 ! vector opt. 360 375 z1_e3t = zfact / fse3t(ji,jj,jk) … … 375 390 IF(lflush) CALL flush(numout) 376 391 ENDIF 377 IF(nn_timing == 2) CALL timing_start('iom_rstput') 392 IF(nn_timing == 2) CALL timing_start('iom_rstput') 378 393 CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b' , qsr_hc ) 379 CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev ) ! default definition in sbcssm 394 CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev ) ! default definition in sbcssm 380 395 IF(nn_timing == 2) CALL timing_stop('iom_rstput') 381 396 ! … … 385 400 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 386 401 CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt ) 387 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt ) 402 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt ) 388 403 ENDIF 389 404 ! ! print mean trends (used for debugging) 390 405 IF(ln_ctl) CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' qsr - Ta: ', mask1=tmask, clinfo3='tra-ta' ) 391 406 ! 392 CALL wrk_dealloc( jpi, jpj, zekb, zekg, zekr ) 393 CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d ) 407 CALL wrk_dealloc( jpi, jpj, zekb, zekg, zekr ) 408 CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d ) 394 409 ! 395 410 IF( nn_timing == 1 ) CALL timing_stop('tra_qsr') … … 407 422 !! from two length scale of penetration (rn_si0,rn_si1) and a ratio 408 423 !! (rn_abs). These parameters are read in the namtra_qsr namelist. The 409 !! default values correspond to clear water (type I in Jerlov' 424 !! default values correspond to clear water (type I in Jerlov' 410 425 !! (1968) classification. 411 426 !! called by tra_qsr at the first timestep (nit000) … … 434 449 IF( nn_timing == 1 ) CALL timing_start('tra_qsr_init') 435 450 ! 436 CALL wrk_alloc( jpi, jpj, zekb, zekg, zekr ) 437 CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 451 CALL wrk_alloc( jpi, jpj, zekb, zekg, zekr ) 452 CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 438 453 ! 439 454 … … 465 480 466 481 IF( ln_traqsr ) THEN ! control consistency 467 ! 482 ! 468 483 IF( .NOT.lk_qsr_bio .AND. ln_qsr_bio ) THEN 469 484 CALL ctl_warn( 'No bio model : force ln_qsr_bio = FALSE ' ) … … 480 495 & ' 2 bands, 3 RGB bands or bio-model light penetration' ) 481 496 ! 482 IF( ln_qsr_rgb .AND. nn_chldta == 0 ) nqsr = 1 497 IF( ln_qsr_rgb .AND. nn_chldta == 0 ) nqsr = 1 483 498 IF( ln_qsr_rgb .AND. nn_chldta == 1 ) nqsr = 2 484 499 IF( ln_qsr_rgb .AND. nn_chldta == 2 ) nqsr = 3 … … 498 513 ENDIF 499 514 ! ! ===================================== ! 500 IF( ln_traqsr ) THEN ! Initialisation of Light Penetration ! 515 IF( ln_traqsr ) THEN ! Initialisation of Light Penetration ! 501 516 ! ! ===================================== ! 502 517 ! 518 ALLOCATE( xsi0r(jpi,jpj) ) 503 519 xsi0r = 1.e0 / rn_si0 504 520 xsi1r = 1.e0 / rn_si1 … … 539 555 zchl = 0.05 ! constant chlorophyll 540 556 irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) 541 zekb(:,:) = rkrgb(1,irgb) ! Separation in R-G-B depending of the chlorophyll 557 zekb(:,:) = rkrgb(1,irgb) ! Separation in R-G-B depending of the chlorophyll 542 558 zekg(:,:) = rkrgb(2,irgb) 543 559 zekr(:,:) = rkrgb(3,irgb) … … 546 562 ze0(:,:,1) = rn_abs 547 563 ze1(:,:,1) = zcoef 548 ze2(:,:,1) = zcoef 564 ze2(:,:,1) = zcoef 549 565 ze3(:,:,1) = zcoef 550 566 zea(:,:,1) = tmask(:,:,1) ! = ( ze0+ze1+z2+ze3 ) * tmask 551 567 552 568 DO jk = 2, nksr+1 553 569 !CDIR NOVERRCHK 554 570 DO jj = 1, jpj 555 !CDIR NOVERRCHK 571 !CDIR NOVERRCHK 556 572 DO ji = 1, jpi 557 zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * xsi0r 573 zc0 = ze0(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * xsi0r(ji,jj) ) 558 574 zc1 = ze1(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekb(ji,jj) ) 559 575 zc2 = ze2(ji,jj,jk-1) * EXP( - e3t_0(ji,jj,jk-1) * zekg(ji,jj) ) … … 566 582 END DO 567 583 END DO 568 END DO 584 END DO 569 585 ! 570 586 DO jk = 1, nksr … … 600 616 DO jj = 1, jpj ! top 400 meters 601 617 DO ji = 1, jpi 602 zc0 = zz0 * EXP( -fsdepw(ji,jj,jk )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk )*xsi1r )603 zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r )604 etot3(ji,jj,jk) = ( zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1) ) * tmask(ji,jj,1) 618 zc0 = zz0 * EXP( -fsdepw(ji,jj,jk )*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk )*xsi1r ) 619 zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r(ji,jj) ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r ) 620 etot3(ji,jj,jk) = ( zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1) ) * tmask(ji,jj,1) 605 621 END DO 606 622 END DO … … 611 627 ENDIF 612 628 ! ! ===================================== ! 613 ELSE ! No light penetration ! 629 ELSE ! No light penetration ! 614 630 ! ! ===================================== ! 615 631 IF(lwp) THEN … … 617 633 WRITE(numout,*) 'tra_qsr_init : NO solar flux penetration' 618 634 WRITE(numout,*) '~~~~~~~~~~~~' 619 IF(lflush) CALL flush(numout) 635 IF(lflush) CALL flush(numout) 620 636 ENDIF 621 637 ENDIF … … 630 646 ENDIF 631 647 ! 632 CALL wrk_dealloc( jpi, jpj, zekb, zekg, zekr ) 633 CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 648 CALL wrk_dealloc( jpi, jpj, zekb, zekg, zekr ) 649 CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 634 650 ! 635 651 IF( nn_timing == 1 ) CALL timing_stop('tra_qsr_init')
Note: See TracChangeset
for help on using the changeset viewer.
