Changeset 8447


Ignore:
Timestamp:
2017-08-18T18:16:42+02:00 (3 years ago)
Author:
anaguiar
Message:

Merging with the dev_r5518_pcbias_ipc@8446 branch

Location:
branches/UKMO/dev_r5518_GO6_package/NEMOGCM
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_GO6_package/NEMOGCM/CONFIG/SHARED/namelist_ref

    r8427 r8447  
    13631363   nn_bias_itwrt  = 15 
    13641364   ln_itdecay     = .FALSE. 
    1365 / 
     1365   ln_incpc       = .FALSE. 
     1366/ 
  • branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/bias.F90

    r8400 r8447  
    7676   USE oce, ONLY: & 
    7777      & tsb, tsn, tsa, & 
    78       & rhop,          & 
     78      & rhop,  & 
    7979      & gtsu, gtsv 
    8080   USE dynhpg, ONLY:   & 
     
    134134      & fb_p_ofl,             &  !: parition of bias in P for ofl bias term 
    135135      & fctamp,               &  !: amplification factor for T if inertial 
    136       & rn_maxlat_bias,         &  !: Max lat for latitudinal ramp 
    137       & rn_minlat_bias             !: Min lat for latitudinal ramp 
     136      & rn_maxlat_bias,       &  !: Max lat for latitudinal ramp 
     137      & rn_minlat_bias,       &  !: Min lat for latitudinal ramp 
     138      & zwgt,                 &  !: weight for IPC 
     139      & ztscale                  !: decay rate for IPC 
    138140 
    139141   LOGICAL,  PRIVATE :: lalloc 
     
    148150      & sbias_rlx_out, &   !: Output salinity bias field 
    149151      & tbias_p_out,   &   !: Output temperature bias field for P correction 
    150       & sbias_p_out        !: Output salinity bias field for P correction 
     152      & sbias_p_out,   &   !: Output salinity bias field for P correction 
     153      & tbias_i_out,   &   !: Output temperature bias field for incremental P correction 
     154      & sbias_i_out,   &   !: Output salinity bias field for incremental P correction 
     155      & tbias_asm_stscale, &   !: Short time scale temperature bias field 
     156      & sbias_asm_stscale, &   !: Short time scale salinity bias field 
     157      & tbias_asm_stscale_out, &   !: Short time scale temperature bias output field 
     158      & sbias_asm_stscale_out   !: Short time scale salinity bias output field 
    151159 
    152160   INTEGER, PRIVATE :: nn_lat_ramp     ! choice of latitude dependent ramp 
     
    154162                   ! 1:inertial ramp, 2:linear ramp, else:no ramp 
    155163   LOGICAL, PRIVATE :: ln_bsyncro      ! syncronous or assincrous bias correction    
    156    LOGICAL, PRIVATE :: ln_itdecay      ! evolve bias correction at every time step.    
    157  
    158    REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: fbcoef   
     164   LOGICAL, PRIVATE :: ln_itdecay      ! evolve bias correction at every time step.   
     165   LOGICAL, PRIVATE :: ln_incpc        ! incremental pressure correction plus pressure correction 
     166   LOGICAL, PRIVATE :: ln_incpc_only        ! incremental pressure correction only 
     167 
     168   REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: fbcoef 
     169   REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: fbcoef_stscale 
    159170 
    160171   INTEGER, PRIVATE ::  & 
     
    200211      !! * Local declarations 
    201212 
    202       CHARACTER(len=100) ::  cn_dir       ! dir for location ofline bias 
     213      CHARACTER(len=100) ::  cn_dir          ! dir for location ofline bias 
    203214      INTEGER            ::  ierror 
    204       INTEGER            ::  ios          ! Local integer output status for namelist read 
    205       REAL(wp)           ::  eft_rlx,  &  ! efolding time (bias memory) in days 
    206          &                   eft_asm,  &  !     "      " 
     215      INTEGER            ::  ios             ! Local integer output status for namelist read 
     216      REAL(wp)           ::  eft_rlx,  &     ! efolding time (bias memory) in days 
     217         &                   eft_asm,  &     !     "      " 
    207218         &                   log2,     & 
    208          &                   lenscl_bias     !lengthscale of the pressure bias decay between minlat and maxlat. 
     219         &                   lenscl_bias, &  ! lengthscale pressure bias decay between minlat and maxlat. 
     220         &                   minlat_bias, &  ! used in ipc 
     221         &                   maxlat_bias     ! used in ipc 
    209222       
    210223      NAMELIST/nambias/ ln_bias, ln_bias_asm, ln_bias_rlx, ln_bias_ofl,   & 
     
    215228         & cn_bias_tot, cn_bias_asm, cn_dir, sn_tbias_ofl, sn_sbias_ofl,  & 
    216229         & ln_bsyncro, fctamp, rn_maxlat_bias, rn_minlat_bias,            & 
    217          & nn_bias_itwrt, ln_itdecay 
     230         & nn_bias_itwrt, ln_itdecay, ln_incpc,ln_incpc_only, ztscale, zwgt 
    218231  
    219232 
     
    274287         WRITE(numout,*) '     time step for writing bias fld nn_bias_itwrt  = ',nn_bias_itwrt 
    275288    WRITE(numout,*) '     evolve pcbias at each timestep ln_itdecay     = ',ln_itdecay 
     289    WRITE(numout,*) '     incremental press. correction  ln_incpc       = ',ln_incpc 
     290    WRITE(numout,*) '     incremental press. correction only  ln_incpc_only       = ',ln_incpc_only 
     291    WRITE(numout,*) '     incremental press. correction  zwgt           = ',zwgt 
     292    WRITE(numout,*) '     incremental press. correction  ztscale        = ',ztscale 
    276293         WRITE(numout,*) '  Parameters for parition of bias term ' 
    277294         WRITE(numout,*) '                                    fb_t_rlx       = ',fb_t_rlx 
     
    296313            &   CALL ctl_stop (' lk_dtatem, lk_dtasal and lk_tradmp need to be true with ln_bias_rlx' ) 
    297314 
     315         IF ( (.NOT. ln_itdecay) .AND. ln_incpc) &    
     316            &   CALL ctl_stop (' if you set ln_incpc to .true. then you need to set ln_itdecay to .true. as well' ) 
     317 
     318         IF ( (.NOT. ln_incpc) .AND. ln_incpc_only) &    
     319            &   CALL ctl_stop (' if you set ln_incpc_only to .true. then you need to set ln_incpc to .true. as well' ) 
     320          
     321         WRITE(numout,*) '     time step is    = ',rdt,'you choose to write pcbias at nn_bias_itwrt  = ',nn_bias_itwrt,'and end of iau is rday/rdt=',rday/rdt  
    298322      ENDIF 
    299323      IF( .NOT. ln_bias ) RETURN 
     
    305329            &      tbias_p(jpi,jpj,jpk), & 
    306330            &      sbias_p(jpi,jpj,jpk), & 
     331            &      tbias_i(jpi,jpj,jpk), & 
     332            &      sbias_i(jpi,jpj,jpk), & 
    307333            &      rhd_pc(jpi,jpj,jpk) , & 
    308334            &      gru_pc(jpi,jpj)     , & 
    309335            &      grv_pc(jpi,jpj)       ) 
    310336 
    311          ALLOCATE( fbcoef(jpi,jpj)      ) 
     337         ALLOCATE( fbcoef(jpi,jpj), fbcoef_stscale(jpi,jpj) ) 
    312338          
    313339         IF( ln_bias_asm ) ALLOCATE(  tbias_asm(jpi,jpj,jpk),     & 
     
    322348                                      tbias_rlx_out(jpi,jpj,jpk), &   
    323349                                      sbias_rlx_out(jpi,jpj,jpk)  ) 
     350 
     351         IF( ln_incpc )    ALLOCATE(  tbias_asm_stscale(jpi,jpj,jpk), & 
     352            &                         sbias_asm_stscale(jpi,jpj,jpk), & 
     353                                      tbias_asm_stscale_out(jpi,jpj,jpk), & 
     354                                      sbias_asm_stscale_out(jpi,jpj,jpk), & 
     355                                      tbias_i_out(jpi,jpj,jpk),   &   
     356                                      sbias_i_out(jpi,jpj,jpk)   ) 
     357 
    324358         lalloc = .TRUE. 
    325359 
     
    360394         tbias_p(:,:,:) = 0.0_wp 
    361395         sbias_p(:,:,:) = 0.0_wp 
     396         tbias_i(:,:,:) = 0.0_wp 
     397         sbias_i(:,:,:) = 0.0_wp 
    362398         gru_pc(:,:)    = 0.0_wp 
    363399         grv_pc(:,:)    = 0.0_wp 
     400          
    364401         IF ( ln_bias_rlx ) THEN 
    365402            tbias_rlx(:,:,:) = 0.0_wp 
    366403            sbias_rlx(:,:,:) = 0.0_wp 
    367404         ENDIF 
     405 
    368406         IF ( ln_bias_asm ) THEN   !now rlx and asm bias in same file 
    369            tbias_asm(:,:,:) = 0.0_wp 
    370            sbias_asm(:,:,:) = 0.0_wp 
    371          ENDIF 
     407            tbias_asm(:,:,:) = 0.0_wp 
     408            sbias_asm(:,:,:) = 0.0_wp 
     409            tbias_asm_out(:,:,:) = 0.0_wp 
     410            sbias_asm_out(:,:,:) = 0.0_wp 
     411         ENDIF 
     412 
     413         IF ( ln_incpc ) THEN   !incr pressure correction 
     414            tbias_asm_stscale(:,:,:) = 0.0_wp 
     415            sbias_asm_stscale(:,:,:) = 0.0_wp 
     416            tbias_asm_stscale_out(:,:,:) = 0.0_wp 
     417            sbias_asm_stscale_out(:,:,:) = 0.0_wp 
     418         ENDIF 
     419 
     420 
    372421         numbias_tot    = 0 
    373422         ! Get bias from file and prevent fail if the file does not exist 
     
    404453               ENDIF 
    405454            ENDIF 
     455 
     456            IF ( ln_incpc .and. .not.ln_bsyncro ) THEN 
     457               IF(lwp) WRITE(numout,*) 'Reading short time scale bias correction fields for bias asm from file ',cn_bias_tot 
     458               IF( iom_varid( numbias_tot, 'tbias_asm_stscale' ) > 0 ) THEN 
     459                  ! Get the T and S bias data 
     460                  CALL iom_get( numbias_tot, jpdom_autoglo, 'tbias_asm_stscale', tbias_asm_stscale ) 
     461                  CALL iom_get( numbias_tot, jpdom_autoglo, 'sbias_asm_stscale', sbias_asm_stscale ) 
     462               ELSE 
     463                  CALL ctl_stop( 'Short time scale bias assim variables not found in ',cn_bias_tot ) 
     464               ENDIF 
     465            ENDIF   
     466 
     467 
    406468 
    407469            ! Close the file 
     
    484546         fbcoef(:,:) = 0.0_wp 
    485547         fctamp      = 0.0_wp 
    486       ENDIF 
     548         fbcoef_stscale(:,:) = 0.0_wp 
     549      ENDIF 
     550 
     551 
     552      IF ( ln_incpc) THEN 
     553         minlat_bias = 3.0_wp 
     554         maxlat_bias = 8.0_wp   
     555         WHERE ( abs( gphit(:,:) ) <= minlat_bias ) 
     556            fbcoef_stscale(:,:)=0._wp 
     557         ELSEWHERE ( abs( gphit(:,:) ) >= maxlat_bias )  
     558            fbcoef_stscale(:,:)=1._wp 
     559         ELSEWHERE 
     560            fbcoef_stscale(:,:)=1._wp - ((maxlat_bias - abs( gphit(:,:)))/(maxlat_bias-minlat_bias)) 
     561         ENDWHERE    
     562      ENDIF  
    487563 
    488564 
     
    517593      REAL(wp)            ::   fb_t_asm_max, fb_t_rlx_max, fb_t_ofl_max 
    518594      REAL(wp)            ::   ztfrac, ztsday 
     595      REAL(wp)            ::   zfrac, zfrac1 ! temporal weights for inst pcbias (names could be changed) 
     596      REAL(wp)            ::   zdecay        ! used in inst pcorr 
    519597      REAL(wp), DIMENSION(jpi,jpj) :: zcof1, zcof2 
    520598 
     
    528606      tbias_p(:,:,:) = 0.0_wp 
    529607      sbias_p(:,:,:) = 0.0_wp 
     608      tbias_i(:,:,:) = 0.0_wp 
     609      sbias_i(:,:,:) = 0.0_wp 
     610 
    530611      IF ( ln_bias_asm ) THEN 
    531612       
     
    535616            &                              fbcoef(:,:)   * fb_t_asm_max ) 
    536617         zcof2(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_p_asm 
    537        
     618 
    538619         IF ( ln_itdecay ) THEN   !decay the pressure correction at each time step 
    539620     
    540        ztsday  = rday / real(rdt) 
     621       ztsday  = rday / REAL(rdt,wp) 
     622 
     623            zdecay = (1-ztscale)**(1/REAL(ztsday,wp)) ! used in ipc 
     624            zfrac1 = zdecay**REAL(kt,wp) ! used in ipc 
     625            IF ( zfrac1 <= 0.0 ) zfrac1 = 0.0_wp 
    541626 
    542627            IF( ln_asmiau .and. ln_trainc ) THEN  ! nudge in increments and decay historical contribution 
    543  
    544628                
    545629               IF ( kt <= nitiaufin ) THEN  ! During IAU calculate the fraction of increments to apply at each time step 
    546630 
    547                   ztfrac = real(kt) / real(nitiaufin)  ! nudging factor during the IAU 
    548  
     631                  ztfrac = REAL(kt,wp) / REAL(nitiaufin,wp)  ! nudging factor during the IAU 
    549632           
    550633                  IF (lwp) THEN 
     
    552635                     WRITE(numout,*) '~~~~~~~~~~~~' 
    553636                     WRITE(numout,* ) "proportion of  increment applied in pcbias ",ztfrac 
    554                      WRITE(numout,* ) "proportion of  historical pcbias applied ",t_asm_mem**(real(kt)/ztsday) 
     637                     WRITE(numout,* ) "proportion of  historical pcbias applied ",t_asm_mem**(REAL(kt,wp)/ztsday) 
    555638                  ENDIF 
    556             
    557                   DO jk = 1, jpkm1 
     639 
     640                  DO jk = 1, jpkm1   
    558641                     tbias(:,:,jk) = tbias(:,:,jk) +                            & 
    559                      &                ( t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk)  +                    & 
     642                     &                ( t_asm_mem**(REAL(kt,wp)/ztsday) * tbias_asm(:,:,jk)  +                    & 
    560643                     &                ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof1(:,:) 
    561644                     sbias(:,:,jk) = sbias(:,:,jk) +                            & 
    562                      &               ( t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk)  +                     & 
     645                     &               ( t_asm_mem**(REAL(kt,wp)/ztsday) * sbias_asm(:,:,jk)  +                     & 
    563646                     &               ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof1(:,:) 
     647 
    564648                     tbias_p(:,:,jk) = tbias_p(:,:,jk) +                        & 
    565                      &               ( t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk)  +                     & 
     649                     &               ( t_asm_mem**(REAL(kt,wp)/ztsday) * tbias_asm(:,:,jk)  +                     & 
    566650                     &               ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof2(:,:) 
    567651                     sbias_p(:,:,jk) = sbias_p(:,:,jk) +                        &   
    568                      &               ( t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk)  +                     & 
     652                     &               ( t_asm_mem**(REAL(kt,wp)/ztsday) * sbias_asm(:,:,jk)  +                     & 
    569653                     &               ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof2(:,:) 
    570654                  ENDDO 
    571655 
     656                  IF (ln_incpc) THEN 
     657 
     658                     IF (lwp) THEN 
     659                       WRITE(numout,*) 'tra_bias : bias weights' 
     660                       WRITE(numout,*) '~~~~~~~~~~~~' 
     661                       WRITE(numout,* ) "IPC - proportion of  increment applied in pcbias ",ztfrac 
     662                       WRITE(numout,* ) "IPC - proportion of  historical pcbias applied ",zfrac1 
     663                     ENDIF 
     664 
     665                     DO jk = 1, jpkm1 
     666 
     667                        tbias_i(:,:,jk) =  ( t_bkginc(:,:,jk) * tmask(:,:,jk)* zwgt * ztfrac * (1.0 - fbcoef_stscale(:,:)) )         & 
     668                        &                + ( tbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) ) 
     669                        sbias_i(:,:,jk) =  ( s_bkginc(:,:,jk) * tmask(:,:,jk)* zwgt * ztfrac * (1.0 - fbcoef_stscale(:,:)) )         & 
     670                        &                + ( sbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) ) 
     671 
     672                     ENDDO 
     673 
     674                  ENDIF 
     675 
    572676                  IF ( .not.ln_bsyncro ) THEN  
     677 
    573678                     IF ( kt == nn_bias_itwrt ) THEN 
     679 
    574680                        DO jk = 1, jpkm1 
    575                            tbias_asm_out(:,:,jk) =  t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk)  +       & 
     681                           tbias_asm_out(:,:,jk) =  t_asm_mem**(REAL(kt,wp)/ztsday) * tbias_asm(:,:,jk)  +       & 
    576682                           &                     ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk) 
    577                            sbias_asm_out(:,:,jk) =  t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk)  +       & 
     683                           sbias_asm_out(:,:,jk) =  t_asm_mem**(REAL(kt,wp)/ztsday) * sbias_asm(:,:,jk)  +       & 
    578684                           &                     ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) 
    579685                         END DO 
     686 
     687                        IF ( ln_incpc) THEN 
     688                           DO jk = 1, jpkm1 
     689                              tbias_asm_stscale_out(:,:,jk) = ( t_bkginc(:,:,jk) * tmask(:,:,jk) *  zwgt * ztfrac ) + ( tbias_asm_stscale(:,:,jk) * zfrac1 ) 
     690                              sbias_asm_stscale_out(:,:,jk) = ( s_bkginc(:,:,jk) * tmask(:,:,jk) *  zwgt * ztfrac ) + ( sbias_asm_stscale(:,:,jk) * zfrac1 ) 
     691                           ENDDO 
     692                        ENDIF 
     693 
    580694                     ENDIF 
     695 
    581696                  ENDIF 
    582697 
     
    584699                  IF ( kt == nitiaufin ) THEN 
    585700                     DO jk = 1, jpkm1 
    586                         tbias_asm(:,:,jk) =  t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk)  +       & 
     701                        tbias_asm(:,:,jk) =  t_asm_mem**(REAL(kt,wp)/ztsday) * tbias_asm(:,:,jk)  +       & 
    587702                        &                     ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk) 
    588                         sbias_asm(:,:,jk) =  t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk)  +       & 
     703                        sbias_asm(:,:,jk) =  t_asm_mem**(REAL(kt,wp)/ztsday) * sbias_asm(:,:,jk)  +       & 
    589704                        &                     ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) 
    590705                     END DO 
     706 
     707                     IF (ln_incpc) THEN 
     708                        DO jk = 1, jpkm1 
     709                           tbias_asm_stscale(:,:,jk) = ( t_bkginc(:,:,jk) * tmask(:,:,jk) * zwgt * ztfrac ) + ( tbias_asm_stscale(:,:,jk) * zfrac1 ) 
     710                           sbias_asm_stscale(:,:,jk) = ( s_bkginc(:,:,jk) * tmask(:,:,jk) * zwgt * ztfrac ) + ( sbias_asm_stscale(:,:,jk) * zfrac1 ) 
     711                        ENDDO 
     712                     ENDIF 
     713      
    591714                  ENDIF 
    592715                
    593716               ELSE ! decay pressure correction from combined historical component and increments after IAU 
    594717 
    595                   ztfrac=t_asm_mem**(real(kt-nitiaufin)/real(nitiaufin))  ! decay from end of IAU 
    596  
     718                  ztfrac=t_asm_mem**(REAL(kt-nitiaufin,wp)/REAL(nitiaufin,wp))  ! decay from end of IAU 
     719                   
    597720                  DO jk = 1, jpkm1 
    598721                     tbias(:,:,jk) = tbias(:,:,jk) +                            & 
     
    606729                  ENDDO 
    607730 
    608                   IF (lwp) THEN 
    609                      WRITE(numout,*) 'tra_bias : bias weights' 
    610                      WRITE(numout,*) '~~~~~~~~~~~~' 
    611                      WRITE(numout,* ) "proportion of increments + historical pcbias applied ",ztfrac 
    612                   ENDIF 
    613  
    614                   IF ( ln_trainc .and. .not.ln_bsyncro ) THEN  
    615                      IF ( kt == nn_bias_itwrt ) THEN 
    616                         DO jk = 1, jpkm1 
    617                            tbias_asm_out(:,:,jk) =  ztfrac * tbias_asm(:,:,jk)  
    618                            sbias_asm_out(:,:,jk) =  ztfrac * sbias_asm(:,:,jk)  
    619                          END DO 
    620                      ENDIF 
    621                   ENDIF 
    622  
    623                ENDIF 
    624  
     731                 IF (ln_incpc) THEN 
     732 
     733                   zfrac  = zdecay**REAL(kt-nitiaufin,wp)  
     734                   IF ( zfrac <= 0.0 ) zfrac = 0.0_wp 
     735    
     736                   DO jk = 1, jpkm1 
     737                      tbias_i(:,:,jk) =  tbias_asm_stscale(:,:,jk) * zfrac * (1.0 - fbcoef_stscale(:,:)) 
     738                      sbias_i(:,:,jk) =  sbias_asm_stscale(:,:,jk) * zfrac * (1.0 - fbcoef_stscale(:,:)) 
     739                   ENDDO 
     740 
     741                   IF (lwp) THEN 
     742                      WRITE(numout,*) 'tra_bias : bias weights' 
     743                      WRITE(numout,*) '~~~~~~~~~~~~' 
     744                      WRITE(numout,* ) "IPC - proportion of increments + historical pcbias applied ",zfrac 
     745                   ENDIF 
     746 
     747                 ENDIF 
     748 
     749                 IF (lwp) THEN 
     750                    WRITE(numout,*) 'tra_bias : bias weights' 
     751                    WRITE(numout,*) '~~~~~~~~~~~~' 
     752                    WRITE(numout,* ) "proportion of increments + historical pcbias applied ",ztfrac 
     753                 ENDIF 
     754 
     755                 IF ( ln_trainc .and. .not.ln_bsyncro ) THEN  
     756                    IF ( kt == nn_bias_itwrt ) THEN 
     757                       DO jk = 1, jpkm1 
     758                          tbias_asm_out(:,:,jk) =  ztfrac * tbias_asm(:,:,jk)  
     759                          sbias_asm_out(:,:,jk) =  ztfrac * sbias_asm(:,:,jk)  
     760                       END DO 
     761 
     762                       IF ( ln_incpc ) THEN 
     763                          IF (lwp) WRITE(numout,*) 'after end of IAU - IPC - saving s/tbias_asm_stscale' 
     764                             DO jk = 1, jpkm1 
     765                                tbias_asm_stscale_out(:,:,jk) =  tbias_asm_stscale(:,:,jk) * zfrac 
     766                                sbias_asm_stscale_out(:,:,jk) =  sbias_asm_stscale(:,:,jk) * zfrac 
     767                             ENDDO 
     768                          ENDIF 
     769                       ENDIF 
     770                    ENDIF 
     771                 ENDIF 
    625772 
    626773            ELSE ! no assimilation increments, simply decay pressure correction (e.g for forecasts) 
     774                 ! this occurs at obsoper step as well ( ln_asmiau is false )  
    627775 
    628776               DO jk = 1, jpkm1 
    629777                  tbias(:,:,jk) = tbias(:,:,jk) +                                                         & 
    630                   &               ( t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk) ) * zcof1(:,:) 
     778                  &               ( t_asm_mem**(REAL(kt,wp)/ztsday) * tbias_asm(:,:,jk) ) * zcof1(:,:) 
    631779                  sbias(:,:,jk) = sbias(:,:,jk) +                                                         & 
    632                   &               ( t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk) ) * zcof1(:,:) 
     780                  &               ( t_asm_mem**(REAL(kt,wp)/ztsday) * sbias_asm(:,:,jk) ) * zcof1(:,:) 
    633781                  tbias_p(:,:,jk) = tbias_p(:,:,jk) +                                                     & 
    634                   &               ( t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk) ) * zcof2(:,:) 
     782                  &               ( t_asm_mem**(REAL(kt,wp)/ztsday) * tbias_asm(:,:,jk) ) * zcof2(:,:) 
    635783                  sbias_p(:,:,jk) = sbias_p(:,:,jk) +                                                     &   
    636                   &               ( t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk) ) * zcof2(:,:) 
     784                  &               ( t_asm_mem**(REAL(kt,wp)/ztsday) * sbias_asm(:,:,jk) ) * zcof2(:,:) 
    637785               ENDDO 
    638786 
     
    640788                  WRITE(numout,*) 'tra_bias : bias weights' 
    641789                  WRITE(numout,*) '~~~~~~~~~~~~' 
    642                   WRITE(numout,* ) "proportion of historical pcbias applied ",t_asm_mem**(real(kt)/ztsday) 
     790                  WRITE(numout,* ) "proportion of historical pcbias applied ",t_asm_mem**(REAL(kt,wp)/ztsday) 
    643791               ENDIF 
    644792 
     793               IF (ln_incpc) THEN 
     794                   IF (lwp) WRITE(numout,*) 'obsoper or forecast mode - IPC - computing tbias_i and sbias_i' 
     795                   DO jk = 1, jpkm1 
     796                      tbias_i(:,:,jk) = tbias_i(:,:,jk) + ( tbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) ) 
     797                      sbias_i(:,:,jk) = sbias_i(:,:,jk) + ( sbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) ) 
     798                   ENDDO 
     799                   tbias_i(:,:,:) =tbias_i(:,:,:)*tmask(:,:,:) 
     800                   sbias_i(:,:,:) =sbias_i(:,:,:)*tmask(:,:,:) 
     801               ENDIF 
    645802            ENDIF 
    646803  
     
    731888 
    732889      IF ( kt == nn_bias_itwrt ) THEN 
    733          tbias_p_out(:,:,:)   = tbias_p(:,:,:) 
    734          sbias_p_out(:,:,:)   = sbias_p(:,:,:) 
     890            tbias_p_out(:,:,:)   = tbias_p(:,:,:) 
     891            sbias_p_out(:,:,:)   = sbias_p(:,:,:) 
     892            IF (ln_incpc) THEN 
     893               tbias_i_out(:,:,:)   = tbias_i(:,:,:) 
     894               sbias_i_out(:,:,:)   = sbias_i(:,:,:)             
     895            ENDIF 
    735896      ENDIF 
    736897 
     
    774935         tsw(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) - tbias_p(:,:,:) 
    775936         tsw(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) - sbias_p(:,:,:) 
     937         IF ( ln_incpc ) THEN 
     938            tsw(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) - tbias_p(:,:,:) - tbias_i(:,:,:) 
     939            tsw(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) - sbias_p(:,:,:) - sbias_i(:,:,:) 
     940            IF ( ln_incpc_only ) THEN 
     941               IF(lwp) WRITE(numout,*) 'ln_incpc_only =', ln_incpc_only, 'tsw updated with IPC only and ln_dynhpg_imp = ',ln_dynhpg_imp 
     942               tsw(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) - tbias_i(:,:,:) 
     943               tsw(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) - sbias_i(:,:,:) 
     944            ENDIF 
     945         ENDIF 
    776946      ELSE 
    777947         tsw(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) - tbias_p(:,:,:) 
    778948         tsw(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) - sbias_p(:,:,:) 
    779       ENDIF 
    780  
    781       CALL eos( tsw, rhd_pc, rhop , fsdept_n(:,:,:)) 
     949         IF ( ln_incpc ) THEN 
     950            tsw(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) - tbias_p(:,:,:) - tbias_i(:,:,:) 
     951            tsw(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) - sbias_p(:,:,:) - sbias_i(:,:,:) 
     952            IF ( ln_incpc_only ) THEN 
     953               IF(lwp) WRITE(numout,*) 'ln_incpc_only =', ln_incpc_only, 'tsw updated with IPC only and ln_dynhpg_imp = ',ln_dynhpg_imp 
     954               tsw(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) - tbias_i(:,:,:) 
     955               tsw(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) - sbias_i(:,:,:) 
     956            ENDIF 
     957         ENDIF 
     958      ENDIF 
     959 
     960      IF(lwp) WRITE(numout,*) 'dyn_bias( kt ) calculating rhd_pc, kt =', kt 
     961      CALL eos( tsw, rhd_pc, rhop, fsdept_n(:,:,:) ) 
    782962       
    783       ! is this needed? 
    784963      CALL lbc_lnk( rhd_pc, 'T', 1.0_wp ) 
     964 
    785965 
    786966      ! Partial steps: now horizontal gradient of t,s,rd  
    787967      ! at the bottom ocean level 
     968 
    788969      IF( ln_zps    )  THEN 
    789970         CALL zps_hde( kt, jpts, tsw, gtsu, gtsv,  & 
     
    8561037         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_p'   , sbias_p_out )          
    8571038      ENDIF 
     1039 
     1040      IF ( ln_incpc ) THEN 
     1041         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_asm_stscale' , tbias_asm_stscale_out )    
     1042         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_asm_stscale' , sbias_asm_stscale_out )  
     1043         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_i'   , tbias_i_out ) 
     1044         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_i'   , sbias_i_out )    
     1045      ENDIF 
    8581046       
    8591047      IF( kt == nitend ) THEN 
  • branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/biaspar.F90

    r8400 r8447  
    2727      & tbias, &       !: Temperature bias field for T correction 
    2828      & tbias_p, &     !:  "           "    "    "   P correction 
     29      & tbias_i, &     !:  "           "    "    "   incremental P correction 
    2930      & sbias, &       !: Salinity bias field    for S correction 
    3031      & sbias_p, &     !:  "           "    "    "   P correction 
     32      & sbias_i, &     !:  "           "    "    "   incremental P correction 
    3133      & rhd_pc         !: Press corrtd density from online to use in dyn_hpg 
    3234 
  • branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r8400 r8447  
    113113 
    114114      ENDIF 
    115       ! 
     115 
    116116      SELECT CASE ( nhpg )      ! Hydrostatic pressure gradient computation 
    117117      CASE (  0 )   ;   CALL hpg_zco    ( kt )      ! z-coordinate 
     
    134134      ! 
    135135      IF ( ln_bias .AND. ln_bias_pc_app )  THEN 
     136         IF(lwp) THEN  
     137         WRITE(numout,*) " ! restore original density" 
     138         ENDIF 
    136139         rhd(:,:,:) = z_rhd_st(:,:,:)     ! restore original density 
    137140         gru(:,:)   = z_gru_st(:,:) 
  • branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90

    r8427 r8447  
    457457      END SELECT 
    458458      ! 
     459      CALL lbc_lnk( prd, 'T', 1.0_wp ) 
     460      ! 
    459461      IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', tab3d_2=prhop, clinfo2=' pot : ', ovlap=1, kdim=jpk ) 
    460462      ! 
Note: See TracChangeset for help on using the changeset viewer.