Changeset 7081


Ignore:
Timestamp:
2016-10-25T09:37:22+02:00 (4 years ago)
Author:
isabella
Message:

Initial changes for ipc

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_pcbias_ipc/NEMOGCM/NEMO/OPA_SRC/ASM/bias.F90

    r6337 r7081  
    133133      & fb_p_ofl,             &  !: parition of bias in P for ofl bias term 
    134134      & fctamp,               &  !: amplification factor for T if inertial 
    135       & rn_maxlat_bias,         &  !: Max lat for latitudinal ramp 
    136       & rn_minlat_bias             !: Min lat for latitudinal ramp 
     135      & rn_maxlat_bias,       &  !: Max lat for latitudinal ramp 
     136      & rn_minlat_bias           !: Min lat for latitudinal ramp 
    137137 
    138138   LOGICAL,  PRIVATE :: lalloc 
     
    153153                   ! 1:inertial ramp, 2:linear ramp, else:no ramp 
    154154   LOGICAL, PRIVATE :: ln_bsyncro      ! syncronous or assincrous bias correction    
    155    LOGICAL, PRIVATE :: ln_itdecay      ! evolve bias correction at every time step.    
    156  
    157    REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: fbcoef   
     155   LOGICAL, PRIVATE :: ln_itdecay      ! evolve bias correction at every time step.   
     156   LOGICAL, PRIVATE :: ln_incpc        ! incremental pressure correction 
     157   LOGICAL, PRIVATE :: lat_ramp_incpc  ! ramp for incremental pressure correction 
     158 
     159   REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: fbcoef 
     160   REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: fbcoef_stscale 
    158161 
    159162   INTEGER, PRIVATE ::  & 
     
    214217         & cn_bias_tot, cn_bias_asm, cn_dir, sn_tbias_ofl, sn_sbias_ofl,  & 
    215218         & ln_bsyncro, fctamp, rn_maxlat_bias, rn_minlat_bias,            & 
    216          & nn_bias_itwrt, ln_itdecay 
     219         & nn_bias_itwrt, ln_itdecay, ln_incpc, lat_ramp_incpc 
    217220  
    218221 
     
    273276         WRITE(numout,*) '     time step for writing bias fld nn_bias_itwrt  = ',nn_bias_itwrt 
    274277    WRITE(numout,*) '     evolve pcbias at each timestep ln_itdecay     = ',ln_itdecay 
     278    WRITE(numout,*) '     incremental press. correction  ln_incpc       = ',ln_incpc 
     279         WRITE(numout,*) '     lat ramp for inc.press.corr    lat_ramp_incpc = ',lat_ramp_incpc 
    275280         WRITE(numout,*) '  Parameters for parition of bias term ' 
    276281         WRITE(numout,*) '                                    fb_t_rlx       = ',fb_t_rlx 
     
    304309            &      tbias_p(jpi,jpj,jpk), & 
    305310            &      sbias_p(jpi,jpj,jpk), & 
     311            &      tbias_i(jpi,jpj,jpk), & 
     312            &      sbias_i(jpi,jpj,jpk), & 
    306313            &      rhd_pc(jpi,jpj,jpk) , & 
    307314            &      gru_pc(jpi,jpj)     , & 
    308315            &      grv_pc(jpi,jpj)       ) 
    309316 
    310          ALLOCATE( fbcoef(jpi,jpj)      ) 
     317         ALLOCATE( fbcoef(jpi,jpj), fbcoef_stscale(jpi,jpj) ) 
    311318          
    312319         IF( ln_bias_asm ) ALLOCATE(  tbias_asm(jpi,jpj,jpk),     & 
     
    321328                                      tbias_rlx_out(jpi,jpj,jpk), &   
    322329                                      sbias_rlx_out(jpi,jpj,jpk)  ) 
     330 
     331         IF( ln_incpc )    ALLOCATE(  tbias_asm_stscale(jpi,jpj,jpk), & 
     332            &                         sbias_asm_stscale(jpi,jpj,jpk)) 
     333 
    323334         lalloc = .TRUE. 
    324335 
     
    369380           sbias_asm(:,:,:) = 0.0_wp 
    370381         ENDIF 
     382 
     383         IF ( ln_incpc ) THEN   !incr pressure correction 
     384           tbias_asm_stscale(:,:,:) = 0.0_wp 
     385           sbias_asm_stscale(:,:,:) = 0.0_wp 
     386         ENDIF 
     387 
     388 
    371389         numbias_tot    = 0 
    372390         ! Get bias from file and prevent fail if the file does not exist 
     
    402420                  CALL ctl_stop( 'Bias assim variables not found in ',cn_bias_tot ) 
    403421               ENDIF 
     422 
     423 
     424               IF ( ln_incpc ) THEN 
     425                  IF(lwp) WRITE(numout,*) 'Reading short time scale bias correction fields for bias asm from file ',cn_bias_tot 
     426                  IF( iom_varid( numbias_tot, 'tbias_asm_stscale' ) > 0 ) THEN 
     427                     ! Get the T and S bias data 
     428                     CALL iom_get( numbias_tot, jpdom_autoglo, 'tbias_asm_stscale', tbias_asm_stscale ) 
     429                     CALL iom_get( numbias_tot, jpdom_autoglo, 'sbias_asm_stscale', sbias_asm_stscale ) 
     430                  ELSE 
     431                     CALL ctl_stop( 'Short time scale bias assim variables not found in ',cn_bias_tot ) 
     432                  ENDIF 
     433               ENDIF   
     434 
    404435            ENDIF 
     436 
     437 
    405438 
    406439            ! Close the file 
     
    483516         fbcoef(:,:) = 0.0_wp 
    484517         fctamp      = 0.0_wp 
    485       ENDIF 
     518         fbcoef_stscale(:,:) = 0.0_wp 
     519      ENDIF 
     520 
     521 
     522      IF ( lat_ramp_incpc) THEN 
     523       ! not sure if this should be a special case of nn_lat_ramp == 2 
     524         minlat_bias = 3.0_wp 
     525         maxlat_bias = 8.0_wp   
     526         WHERE ( abs( gphit(:,:) ) <= minlat_bias ) 
     527            fbcoef_stscale(:,:)=0._wp 
     528         ELSEWHERE ( abs( gphit(:,:) ) >= maxlat_bias )  
     529            fbcoef_stscale(:,:)=1._wp 
     530         ELSEWHERE 
     531            fbcoef_stscale(:,:)=1._wp - ((maxlat_bias - abs( gphit(:,:)))/(maxlat_bias-minlat_bias)) 
     532         ENDWHERE    
     533      ENDIF  
    486534 
    487535 
     
    516564      REAL(wp)            ::   fb_t_asm_max, fb_t_rlx_max, fb_t_ofl_max 
    517565      REAL(wp)            ::   ztfrac, ztsday 
     566      REAL(wp)            ::   zfrac, zfrac1 ! temporal weights for inst pcbias (names could be changed) 
     567      REAL(wp)            ::   ztscale       ! reduce the inst pcbias by this amount per 24 hours 
     568      REAL(wp)            ::   zwgt          ! Weight for the inst pcorr term 
     569      REAL(wp)            ::   zdecay        ! used in inst pcorr 
    518570      REAL(wp), DIMENSION(jpi,jpj) :: zcof1, zcof2 
    519571 
     
    527579      tbias_p(:,:,:) = 0.0_wp 
    528580      sbias_p(:,:,:) = 0.0_wp 
     581 
     582      ztscale = 0.1_wp 
     583      zwgt    = 1.0_wp 
     584      tbias_i(:,:,:) = 0.0_wp 
     585      sbias_i(:,:,:) = 0.0_wp 
     586 
    529587      IF ( ln_bias_asm ) THEN 
    530588       
     
    539597       ztsday  = rday / real(rdt) 
    540598 
     599            zdecay = (1-ztscale)**(1/real(ztsday)) ! used in ipc 
     600            zfrac1 = max(0.0_wp, zdecay**real(kt)) ! used in ipc 
     601 
     602 
    541603            IF( ln_asmiau .and. ln_trainc ) THEN  ! nudge in increments and decay historical contribution 
    542604 
     
    545607 
    546608                  ztfrac = real(kt) / real(nitiaufin)  ! nudging factor during the IAU 
    547  
     609                   
    548610           
    549611                  IF (lwp) THEN 
     
    569631                  ENDDO 
    570632 
     633 
     634                  IF (ln_incpc) THEN 
     635 
     636                     DO jk = 1, jpkm1 
     637 
     638                        tbias_i(:,:,jk) = tbias_i(:,:,jk) +                         &             ! not sure if it should sum to previous               
     639                        &               ( t_bkginc(:,:,jk) * zwgt * ztfrac * (1.0 - fbcoef_stscale(:,:)) )         & 
     640                        &                - ( tbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) ) 
     641                        sbias_i(:,:,jk) = sbias_i(:,:,jk) +                         &                             
     642                        &               ( s_bkginc(:,:,jk) * zwgt * ztfrac * (1.0 - fbcoef_stscale(:,:)) )         & 
     643                        &                - ( sbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) ) 
     644 
     645                     ENDDO 
     646 
     647                  ENDIF 
     648 
    571649                  IF ( .not.ln_bsyncro ) THEN  
    572650                     IF ( kt == nn_bias_itwrt ) THEN 
     
    593671 
    594672                  ztfrac=t_asm_mem**(real(kt-nitiaufin)/real(nitiaufin))  ! decay from end of IAU 
    595  
     673                   
     674                   
    596675                  DO jk = 1, jpkm1 
    597676                     tbias(:,:,jk) = tbias(:,:,jk) +                            & 
     
    603682                     sbias_p(:,:,jk) = sbias_p(:,:,jk) +                        &   
    604683                     &               ( ztfrac * sbias_asm(:,:,jk) ) * zcof2(:,:) 
     684 
    605685                  ENDDO 
     686 
     687                 IF (ln_incpc) THEN 
     688 
     689                   zfrac  = max(0.0_wp, zdecay**real(kt-nitiaufin) ) 
     690                    
     691                   DO jk = 1, jpkm1 
     692                      tbias_i(:,:,jk) = tbias_i(:,:,jk) +                         &                ! not sure if it should sum to previous           
     693                      &               ( t_bkginc(:,:,jk) * zwgt * zfrac * (1.0 - fbcoef_stscale(:,:)) )         & 
     694                      &                - ( tbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) ) 
     695                      sbias_i(:,:,jk) = sbias_i(:,:,jk) +                         &                            ! not sure if it should sum to previous 
     696                      &               ( s_bkginc(:,:,jk) * zwgt * zfrac * (1.0 - fbcoef_stscale(:,:)) )         & 
     697                      &                - ( sbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) ) 
     698                    ENDDO 
     699 
     700                 ENDIF 
    606701 
    607702                  IF (lwp) THEN 
     
    773868         tsw(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) - tbias_p(:,:,:) 
    774869         tsw(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) - sbias_p(:,:,:) 
     870 
     871         IF (ln_incpc) THEN 
     872            tsw(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) -tbias_i(:,:,:) - tbias_p(:,:,:) 
     873            tsw(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) -sbias_i(:,:,:) - sbias_p(:,:,:) 
     874         ENDIF 
     875 
    775876      ELSE 
    776877         tsw(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) - tbias_p(:,:,:) 
    777878         tsw(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) - sbias_p(:,:,:) 
     879 
     880         IF (ln_incpc) THEN 
     881            tsw(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) - tbias_i(:,:,:) - tbias_p(:,:,:) 
     882            tsw(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) - sbias_i(:,:,:) - sbias_p(:,:,:) 
     883         ENDIF 
     884 
    778885      ENDIF 
    779886 
Note: See TracChangeset for help on using the changeset viewer.