Changeset 6300 for branches/UKMO/dev_r5518_pcbias
- Timestamp:
- 2016-02-09T12:30:24+01:00 (8 years ago)
- Location:
- branches/UKMO/dev_r5518_pcbias/NEMOGCM
- Files:
-
- 2 edited
-
CONFIG/SHARED/namelist_ref (modified) (1 diff)
-
NEMO/OPA_SRC/ASM/bias.F90 (modified) (14 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_pcbias/NEMOGCM/CONFIG/SHARED/namelist_ref
r6278 r6300 1311 1311 cn_bias_tot = 'pcbias.nc' 1312 1312 bias_time_unit_asm = 86400.0 1313 ln_inertial = .false.1314 / 1313 nn_lat_ramp = 0 1314 / -
branches/UKMO/dev_r5518_pcbias/NEMOGCM/NEMO/OPA_SRC/ASM/bias.F90
r6279 r6300 137 137 & fb_p_asm, & !: parition of bias in P for asm bias term 138 138 & fb_p_ofl, & !: parition of bias in P for ofl bias term 139 & fctamp !: amplification factor for T if inertial 139 & fctamp, & !: amplification factor for T if inertial 140 & rn_maxlat, & !: Max lat for latitudinal ramp 141 & rn_minlat !: Min lat for latitudinal ramp 140 142 141 143 LOGICAL, PRIVATE :: lalloc … … 144 146 & sbias_asm, & !: Salinity bias field 145 147 & tbias_rlx, & !: Temperature bias field 146 & sbias_rlx !: Salinity bias field 147 148 LOGICAL, PRIVATE :: ln_inertial ! spatial partion of bias into 149 !adiabatic/diabatic corrections 150 LOGICAL, PRIVATE :: ln_bsyncro ! syncronous or assincrous bias correction 151 148 & sbias_rlx, & !: Salinity bias field 149 & tbias_asm_out, & !: Output temperature bias field 150 & sbias_asm_out, & !: Output salinity bias field 151 & tbias_rlx_out, & !: Output temperature bias field 152 & sbias_rlx_out, & !: Output salinity bias field 153 & tbias_p_out, & !: Output temperature bias field for P correction 154 & sbias_p_out !: Output salinity bias field for P correction 155 156 INTEGER, PRIVATE :: nn_lat_ramp ! choice of latitude dependent ramp 157 ! for the pressure correction. 158 ! 1:inertial ramp, 2:linear ramp, else:no ramp 159 LOGICAL, PRIVATE :: ln_bsyncro ! syncronous or assincrous bias correction 160 LOGICAL, PRIVATE :: ln_itdecay ! evolve bias correction at every time step. 152 161 153 162 REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: fbcoef 154 163 155 164 INTEGER, PRIVATE :: & 156 & numbias_asm, & ! netcdf id of bias file from assim 157 & numbias_tot ! netcdf id of bias file with total bias 165 & numbias_asm, & ! netcdf id of bias file from assim 166 & numbias_tot, & ! netcdf id of bias file with total bias 167 & nn_bias_itwrt ! time step for outputting bias pressure corr 158 168 159 169 CHARACTER(LEN=128), PRIVATE :: & … … 200 210 & eft_asm, & ! " " 201 211 & log2, & 202 & minlat_bias, & !latitude, poleward of which the pressure bias begins to decay.203 & maxlat_bias, & !latitude, poleward of which the pressure bias is zero.204 212 & lenscl_bias !lengthscale of the pressure bias decay between minlat and maxlat. 205 213 … … 207 215 & ln_bias_ts_app, ln_bias_pc_app, & 208 216 & fb_t_asm, fb_t_rlx, fb_t_ofl, fb_p_asm, fb_p_rlx, fb_p_ofl, & 209 & eft_rlx, t_rlx_upd, eft_asm, t_asm_upd, ln_inertial, &217 & eft_rlx, t_rlx_upd, eft_asm, t_asm_upd, nn_lat_ramp, & 210 218 & bias_time_unit_asm, bias_time_unit_rlx, bias_time_unit_ofl, & 211 219 & cn_bias_tot, cn_bias_asm, cn_dir, sn_tbias_ofl, sn_sbias_ofl, & 212 & ln_bsyncro, fctamp 220 & ln_bsyncro, fctamp, rn_maxlat, rn_minlat, nn_bias_itwrt, ln_itdecay 213 221 214 222 … … 244 252 fb_p_asm = 1.0_wp 245 253 fb_p_ofl = 0.0_wp 246 ln_inertial = .FALSE.254 nn_lat_ramp = 0 247 255 ln_bsyncro = .FALSE. 256 ln_itdecay = .FALSE. 257 rn_maxlat = 23.0_wp 258 rn_minlat = 10.0_wp 259 260 IF ( lk_asminc ) THEN 261 nn_bias_itwrt = nitiaufin 262 ELSE 263 nn_bias_itwrt = nitend 264 ENDIF 248 265 249 266 … … 291 308 WRITE(numout,*) ' bias pressure correctn apply ln_bias_pc_app = ',ln_bias_pc_app 292 309 WRITE(numout,*) ' bias pressure correctn apply ln_bias_pc_app = ',ln_bias_pc_app 293 WRITE(numout,*) ' bias applied according to lat ln_inertial = ',ln_inertial 310 WRITE(numout,*) ' lat ramp for bias correction nn_lat_ramp = ',nn_lat_ramp 311 WRITE(numout,*) ' time step for writing bias fld nn_bias_itwrt = ',nn_bias_itwrt 312 WRITE(numout,*) ' evolve pcbias at each timestep ln_itdecay = ',ln_itdecay 294 313 WRITE(numout,*) ' Parameters for parition of bias term ' 295 314 WRITE(numout,*) ' fb_t_rlx = ',fb_t_rlx … … 329 348 ALLOCATE( fbcoef(jpi,jpj) ) 330 349 331 IF( ln_bias_asm ) ALLOCATE( tbias_asm(jpi,jpj,jpk), & 332 & sbias_asm(jpi,jpj,jpk) ) 333 334 IF( ln_bias_rlx ) ALLOCATE( tbias_rlx(jpi,jpj,jpk), & 335 & sbias_rlx(jpi,jpj,jpk) ) 336 350 IF( ln_bias_asm ) ALLOCATE( tbias_asm(jpi,jpj,jpk), & 351 & sbias_asm(jpi,jpj,jpk), & 352 tbias_asm_out(jpi,jpj,jpk), & 353 sbias_asm_out(jpi,jpj,jpk), & 354 tbias_p_out(jpi,jpj,jpk), & 355 sbias_p_out(jpi,jpj,jpk) ) 356 357 IF( ln_bias_rlx ) ALLOCATE( tbias_rlx(jpi,jpj,jpk), & 358 & sbias_rlx(jpi,jpj,jpk), & 359 tbias_rlx_out(jpi,jpj,jpk), & 360 sbias_rlx_out(jpi,jpj,jpk) ) 337 361 lalloc = .TRUE. 338 362 … … 473 497 474 498 !latitudinal dependence of partition coeficients. Adhoc 475 IF ( ln_inertial) THEN499 IF ( nn_lat_ramp == 1 ) THEN 476 500 !!!fbcoef(:,:) = SIN( rad * gphit(:,:) )**2 477 501 !!! Introduce the (also adhoc) FOAM parameterisation for latitudinal dependence. 478 502 !!! This should be added as a different namelist option later. MM. 08/2011. 479 minlat_bias = 10._wp 480 maxlat_bias = 23._wp 481 lenscl_bias = ( maxlat_bias - minlat_bias )*2._wp 482 WHERE ( abs( gphit(:,:) ) <= minlat_bias ) 503 lenscl_bias = ( rn_maxlat - rn_minlat )*2._wp 504 WHERE ( abs( gphit(:,:) ) <= rn_minlat ) 483 505 fbcoef(:,:) = 0._wp 484 ELSEWHERE ( abs( gphit(:,:) ) >= maxlat_bias)506 ELSEWHERE ( abs( gphit(:,:) ) >= rn_maxlat ) 485 507 fbcoef(:,:) = 1._wp 486 508 ELSEWHERE 487 fbcoef(:,:) = 1._wp - exp( -( abs( gphit(:,:) ) - minlat_bias ) & 488 * ( abs( gphit(:,:) ) - minlat_bias ) / lenscl_bias ) 489 ENDWHERE 509 fbcoef(:,:) = 1._wp - exp( -( abs( gphit(:,:) ) - rn_minlat ) & 510 * ( abs( gphit(:,:) ) - rn_minlat ) / lenscl_bias ) 511 ENDWHERE 512 ELSEIF ( nn_lat_ramp == 2 ) THEN 513 ! Use a linear ramp consist with the geostrophic velocity balance ramp in NEMOVAR 514 515 WHERE ( abs( gphit(:,:) ) <= rn_minlat ) 516 fbcoef(:,:) = 0._wp 517 ELSEWHERE ( abs( gphit(:,:) ) >= rn_maxlat ) 518 fbcoef(:,:) = 1._wp 519 ELSEWHERE 520 fbcoef(:,:) = 1._wp - ((rn_maxlat - abs( gphit(:,:)))/(rn_maxlat - rn_minlat)) 521 ENDWHERE 490 522 ELSE 491 523 fbcoef(:,:) = 0.0_wp 492 524 fctamp = 0.0_wp 493 525 ENDIF 526 494 527 495 528 END SUBROUTINE bias_init … … 522 555 REAL(wp) :: tsclf ! time saling factor 523 556 REAL(wp) :: fb_t_asm_max, fb_t_rlx_max, fb_t_ofl_max 557 REAL(wp) :: ztfrac, ztsday 524 558 REAL(wp), DIMENSION(jpi,jpj) :: zcof1, zcof2 525 559 … … 534 568 sbias_p(:,:,:) = 0.0_wp 535 569 IF ( ln_bias_asm ) THEN 570 536 571 tsclf = 1 537 572 IF ( .NOT.ln_bsyncro ) tsclf = rdt / bias_time_unit_asm … … 539 574 & fbcoef(:,:) * fb_t_asm_max ) 540 575 zcof2(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_p_asm 541 DO jk = 1, jpkm1 542 tbias(:,:,jk) = tbias(:,:,jk) + tbias_asm(:,:,jk) * zcof1(:,:) 543 sbias(:,:,jk) = sbias(:,:,jk) + sbias_asm(:,:,jk) * zcof1(:,:) 544 tbias_p(:,:,jk) = tbias_p(:,:,jk) + tbias_asm(:,:,jk) * zcof2(:,:) 545 sbias_p(:,:,jk) = sbias_p(:,:,jk) + sbias_asm(:,:,jk) * zcof2(:,:) 546 END DO 547 IF( lk_asminc .and. ln_trainc .and. .not.ln_bsyncro ) THEN 548 ! if last outer loop (lk_asminc=true and ln_trainc=true). t/sbias_asm 549 ! is updated, only once (end of run) taking into account units. 550 IF ( (kt == nitend) .and. lrst_bias) THEN 551 IF(lwp) WRITE(numout,*)' estimating asm bias at last timestep: ',kt 552 DO jk = 1, jpkm1 553 tbias_asm(:,:,jk) = t_asm_mem * tbias_asm(:,:,jk) + & 554 & t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk) 555 sbias_asm(:,:,jk) = t_asm_mem * sbias_asm(:,:,jk) + & 556 & t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) 557 END DO 576 577 IF ( ln_itdecay ) THEN !decay the pressure correction at each time step 578 579 ztsday = rday / real(rdt) 580 581 IF( lk_asminc ) THEN ! nudge in increments and decay historical contribution 582 583 584 IF ( kt <= nitiaufin ) THEN ! During IAU calculate the fraction of increments to apply at each time step 585 586 ztfrac = real(kt) / real(nitiaufin) ! nudging factor during the IAU 587 588 589 IF (lwp) THEN 590 WRITE(numout,*) 'tra_bias : bias weights' 591 WRITE(numout,*) '~~~~~~~~~~~~' 592 WRITE(numout,* ) "proportion of increment applied in pcbias ",ztfrac 593 WRITE(numout,* ) "proportion of historical pcbias applied ",t_asm_mem**(real(kt)/ztsday) 594 ENDIF 595 596 DO jk = 1, jpkm1 597 tbias(:,:,jk) = tbias(:,:,jk) + & 598 & ( t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk) + & 599 & ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof1(:,:) 600 sbias(:,:,jk) = sbias(:,:,jk) + & 601 & ( t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk) + & 602 & ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof1(:,:) 603 tbias_p(:,:,jk) = tbias_p(:,:,jk) + & 604 & ( t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk) + & 605 & ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof2(:,:) 606 sbias_p(:,:,jk) = sbias_p(:,:,jk) + & 607 & ( t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk) + & 608 & ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof2(:,:) 609 ENDDO 610 611 IF ( ln_trainc .and. .not.ln_bsyncro ) THEN 612 IF ( kt == nn_bias_itwrt ) THEN 613 DO jk = 1, jpkm1 614 tbias_asm_out(:,:,jk) = t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk) + & 615 & ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk) 616 sbias_asm_out(:,:,jk) = t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk) + & 617 & ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) 618 END DO 619 ENDIF 620 ENDIF 621 622 ! update the historical component with the increments at the end of the IAU 623 IF ( kt == nitiaufin ) THEN 624 DO jk = 1, jpkm1 625 tbias_asm(:,:,jk) = t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk) + & 626 & ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk) 627 sbias_asm(:,:,jk) = t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk) + & 628 & ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) 629 END DO 630 ENDIF 631 632 ELSE ! decay pressure correction from combined historical component and increments after IAU 633 634 ztfrac=t_asm_mem**(real(kt-nitiaufin)/real(nitiaufin)) ! decay from end of IAU 635 636 DO jk = 1, jpkm1 637 tbias(:,:,jk) = tbias(:,:,jk) + & 638 & ( ztfrac * tbias_asm(:,:,jk) ) * zcof1(:,:) 639 sbias(:,:,jk) = sbias(:,:,jk) + & 640 & ( ztfrac * sbias_asm(:,:,jk) ) * zcof1(:,:) 641 tbias_p(:,:,jk) = tbias_p(:,:,jk) + & 642 & ( ztfrac * tbias_asm(:,:,jk) ) * zcof2(:,:) 643 sbias_p(:,:,jk) = sbias_p(:,:,jk) + & 644 & ( ztfrac * sbias_asm(:,:,jk) ) * zcof2(:,:) 645 ENDDO 646 647 IF (lwp) THEN 648 WRITE(numout,*) 'tra_bias : bias weights' 649 WRITE(numout,*) '~~~~~~~~~~~~' 650 WRITE(numout,* ) "proportion of increments + historical pcbias applied ",ztfrac 651 ENDIF 652 653 IF ( ln_trainc .and. .not.ln_bsyncro ) THEN 654 IF ( kt == nn_bias_itwrt ) THEN 655 DO jk = 1, jpkm1 656 tbias_asm_out(:,:,jk) = ztfrac * tbias_asm(:,:,jk) 657 sbias_asm_out(:,:,jk) = ztfrac * sbias_asm(:,:,jk) 658 END DO 659 ENDIF 660 ENDIF 661 662 ENDIF 663 664 665 ELSE ! no assimilation increments, simply decay pressure correction (e.g for forecasts) 666 667 DO jk = 1, jpkm1 668 tbias(:,:,jk) = tbias(:,:,jk) + & 669 & ( t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk) ) * zcof1(:,:) 670 sbias(:,:,jk) = sbias(:,:,jk) + & 671 & ( t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk) ) * zcof1(:,:) 672 tbias_p(:,:,jk) = tbias_p(:,:,jk) + & 673 & ( t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk) ) * zcof2(:,:) 674 sbias_p(:,:,jk) = sbias_p(:,:,jk) + & 675 & ( t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk) ) * zcof2(:,:) 676 ENDDO 677 678 IF (lwp) THEN 679 WRITE(numout,*) 'tra_bias : bias weights' 680 WRITE(numout,*) '~~~~~~~~~~~~' 681 WRITE(numout,* ) "proportion of historical pcbias applied ",t_asm_mem**(real(kt)/ztsday) 682 ENDIF 683 558 684 ENDIF 559 ENDIF 685 686 ELSE ! maintain a constant pressure correction 687 688 DO jk = 1, jpkm1 689 tbias(:,:,jk) = tbias(:,:,jk) + tbias_asm(:,:,jk) * zcof1(:,:) 690 sbias(:,:,jk) = sbias(:,:,jk) + sbias_asm(:,:,jk) * zcof1(:,:) 691 tbias_p(:,:,jk) = tbias_p(:,:,jk) + tbias_asm(:,:,jk) * zcof2(:,:) 692 sbias_p(:,:,jk) = sbias_p(:,:,jk) + sbias_asm(:,:,jk) * zcof2(:,:) 693 END DO 694 695 IF( lk_asminc .and. ln_trainc .and. .not.ln_bsyncro ) THEN 696 ! if last outer loop (lk_asminc=true and ln_trainc=true). t/sbias_asm 697 ! is updated, only once (end of run) taking into account units. 698 IF ( (kt == nn_bias_itwrt) .and. lrst_bias) THEN 699 IF(lwp) WRITE(numout,*)' estimating asm bias at last timestep: ',kt 700 DO jk = 1, jpkm1 701 tbias_asm_out(:,:,jk) = t_asm_mem * tbias_asm(:,:,jk) + & 702 & t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk) 703 sbias_asm_out(:,:,jk) = t_asm_mem * sbias_asm(:,:,jk) + & 704 & t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) 705 END DO 706 ENDIF 707 ENDIF 708 709 ENDIF 710 560 711 ENDIF 561 712 … … 577 728 sbias_p(:,:,jk) = sbias_p(:,:,jk) + sbias_rlx(:,:,jk) * zcof2(:,:) 578 729 ENDDO 730 IF ( kt == nn_bias_itwrt ) THEN 731 tbias_rlx_out(:,:,:) = tbias_rlx(:,:,:) 732 sbias_rlx_out(:,:,:) = sbias_rlx(:,:,:) 733 ENDIF 579 734 ENDIF 580 735 #endif … … 612 767 CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1.0_wp ) 613 768 769 ENDIF 770 771 IF ( kt == nn_bias_itwrt ) THEN 772 tbias_p_out(:,:,:) = tbias_p(:,:,:) 773 sbias_p_out(:,:,:) = sbias_p(:,:,:) 614 774 ENDIF 615 775 … … 725 885 726 886 IF ( ln_bias_rlx ) THEN 727 CALL iom_rstput( kt, nitend, numbias_tot, 'tbias_rlx' , tbias_rlx)728 CALL iom_rstput( kt, nitend, numbias_tot, 'sbias_rlx' , sbias_rlx)887 CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_rlx' , tbias_rlx_out ) 888 CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_rlx' , sbias_rlx_out ) 729 889 ENDIF 730 890 731 891 IF ( ln_bias_asm ) THEN 732 CALL iom_rstput( kt, nitend, numbias_tot, 'tbias_asm' , tbias_asm)733 CALL iom_rstput( kt, nitend, numbias_tot, 'sbias_asm' , sbias_asm)734 CALL iom_rstput( kt, nitend, numbias_tot, 'tbias_p' , tbias_p)735 CALL iom_rstput( kt, nitend, numbias_tot, 'sbias_p' , sbias_p)892 CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_asm' , tbias_asm_out ) 893 CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_asm' , sbias_asm_out ) 894 CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_p' , tbias_p_out ) 895 CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_p' , sbias_p_out ) 736 896 ENDIF 737 897
Note: See TracChangeset
for help on using the changeset viewer.
