Changeset 11478


Ignore:
Timestamp:
2019-08-28T11:30:41+02:00 (22 months ago)
Author:
mattmartin
Message:

Include initial version of updated pressure correction by Mike Bell.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14_updated_pcbias/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r8447 r11478  
    1717   !!                 !           (A. Coward) suppression of hel, wdj and rot options 
    1818   !!            3.6  !  2014-11  (P. Mathiot) hpg_isf: original code for ice shelf cavity 
     19   !!            3.6  !  2019-08  (M. Bell) Revisions to bias pressure correction  
    1920   !!---------------------------------------------------------------------- 
    2021 
     
    8485      !!---------------------------------------------------------------------- 
    8586      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    86       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv 
    87       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z_rhd_st  ! tmp density storage for pressure corr 
    88       REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   z_gru_st  ! tmp ua trends storage for pressure corr 
    89       REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   z_grv_st  ! tmp va trends storage for pressure corr 
     87      INTEGER                                   ::   ji, jj, jk   ! dummy loop indices 
     88      INTEGER                                   ::   iku, ikv          ! k indices for bottom level at u and v points 
     89      REAL(wp), POINTER, DIMENSION(:,:,:)       ::   ztrdu, ztrdv 
     90      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z_rhd_st   ! tmp density storage for pressure corr 
     91      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z_ua, z_va ! tmp store for ua and va including hpg but not pressure correction  
     92      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z_ua_bpc, z_va_bpc ! ua calculated with bias pressure correction   
     93 
     94      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   z_gru_st, z_grv_st         ! tmp ua and va trends storage for pressure corr 
     95      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   z_ua_bpc_bot, z_va_bpc_bot ! bias pc fields calculated at the ocean bottom 
     96       
    9097      !!---------------------------------------------------------------------- 
    9198      ! 
     
    98105      ENDIF 
    99106      ! 
     107 
    100108      IF ( ln_bias .AND. ln_bias_pc_app ) THEN 
    101109 
     
    103111         ALLOCATE( z_rhd_st(jpi,jpj,jpk), & 
    104112            &      z_gru_st(jpi,jpj),     & 
    105             &      z_grv_st(jpi,jpj)      ) 
    106  
    107          z_rhd_st(:,:,:) = rhd(:,:,:)     ! store orig density  
    108          rhd(:,:,:)      = rhd_pc(:,:,:)  ! use pressure corrected density 
    109          z_gru_st(:,:)   = gru(:,:) 
    110          gru(:,:)        = gru_pc(:,:) 
    111          z_grv_st(:,:)   = grv(:,:) 
    112          grv(:,:)        = grv_pc(:,:) 
    113  
    114       ENDIF 
     113            &      z_grv_st(jpi,jpj),     & 
     114       &      z_ua(jpi,jpj,jpk),     z_va(jpi,jpj,jpk),     &   
     115       &      z_ua_bpc(jpi,jpj,jpk), z_va_bpc(jpi,jpj,jpk), &  
     116            &      z_ua_bpc_bot(jpi,jpj), z_va_bpc_bot(jpi,jpj)  &  
     117       &            ) 
     118 
     119! save the original acceleration trends (z_ua_bpc, z_va_bpc are used as temporary storage) 
     120         z_ua_bpc(:,:,:)     = ua(:,:,:) 
     121    z_va_bpc(:,:,:)     = va(:,:,:) 
     122      END IF  
    115123 
    116124      SELECT CASE ( nhpg )      ! Hydrostatic pressure gradient computation 
     
    123131      END SELECT 
    124132      ! 
     133 
     134      IF ( ln_bias .AND. ln_bias_pc_app ) THEN 
     135 
     136         z_rhd_st(:,:,:) = rhd(:,:,:)     ! store orig density  
     137         rhd(:,:,:)      = rhd_pc(:,:,:)  ! use pressure corrected density 
     138         z_gru_st(:,:)   = gru(:,:) 
     139         gru(:,:)        = gru_pc(:,:) 
     140         z_grv_st(:,:)   = grv(:,:) 
     141         grv(:,:)        = grv_pc(:,:) 
     142 
     143         ! save the acceleration trends including hpg field but calculated without the bpc fields 
     144         z_ua(:,:,:)     = ua(:,:,:) 
     145    z_va(:,:,:)     = va(:,:,:) 
     146 
     147         ! reset the acceleration trends to their original values  
     148         ua(:,:,:)       = z_ua_bpc(:,:,:)     
     149    va(:,:,:)       = z_va_bpc(:,:,:)   
     150 
     151         ! re-calculate the horizontal pressure gradients with the bpc fields  
     152         SELECT CASE ( nhpg )      ! Hydrostatic pressure gradient computation 
     153         CASE (  0 )   ;   CALL hpg_zco    ( kt )      ! z-coordinate 
     154         CASE (  1 )   ;   CALL hpg_zps    ( kt )      ! z-coordinate plus partial steps (interpolation) 
     155         CASE (  2 )   ;   CALL hpg_sco    ( kt )      ! s-coordinate (standard jacobian formulation) 
     156         CASE (  3 )   ;   CALL hpg_djc    ( kt )      ! s-coordinate (Density Jacobian with Cubic polynomial) 
     157         CASE (  4 )   ;   CALL hpg_prj    ( kt )      ! s-coordinate (Pressure Jacobian scheme) 
     158         CASE (  5 )   ;   CALL hpg_isf    ( kt )      ! s-coordinate similar to sco modify for ice shelf 
     159         END SELECT 
     160 
     161         ! calculate the bpc contribution to ua and va 
     162         z_ua_bpc(:,:,:) = ua(:,:,:) - z_ua(:,:,:) 
     163         z_va_bpc(:,:,:) = va(:,:,:) - z_va(:,:,:) 
     164 
     165         ! calculate the bpc contribution to ua and va at the bottom  
     166         DO jj = 2, jpjm1 
     167            DO ji = 2, jpim1 
     168               iku = mbku(ji,jj) 
     169               ikv = mbkv(ji,jj)  
     170               z_ua_bpc_bot(ji,jj) = z_ua_bpc(ji,jj,iku) 
     171               z_va_bpc_bot(ji,jj) = z_va_bpc(ji,jj,ikv) 
     172            END DO  ! ji 
     173         END DO ! jj 
     174 
     175         ! subtract off the bottom values of bpc contribution to ua and va  
     176    DO jk = 1, jpk - 1 
     177           z_ua_bpc(:,:,jk) = z_ua_bpc(:,:,jk) - z_ua_bpc_bot(:,:) 
     178           z_va_bpc(:,:,jk) = z_va_bpc(:,:,jk) - z_va_bpc_bot(:,:) 
     179    END DO  
     180 
     181         ! calculate ua using the original hpg (z_ua) and the bias hpg with the bottom pressure gradients subtracted off  
     182         ua(:,:,:) = z_ua(:,:,:) + z_ua_bpc(:,:,:) 
     183         va(:,:,:) = z_va(:,:,:) + z_va_bpc(:,:,:) 
     184 
     185         IF(lwp) THEN  
     186           WRITE(numout,*) " ! restore original density" 
     187         ENDIF 
     188 
     189         ! restore original density, gru and grv fields  
     190         rhd(:,:,:) = z_rhd_st(:,:,:)      
     191         gru(:,:)   = z_gru_st(:,:) 
     192         grv(:,:)   = z_grv_st(:,:) 
     193 
     194         !Deallocate tempory variables 
     195         DEALLOCATE( z_rhd_st,   z_gru_st, z_grv_st,   & 
     196            &        z_ua, z_va, z_ua_bpc, z_va_bpc,   & 
     197       &        z_ua_bpc_bot, z_va_bpc_bot        & 
     198            &                                          ) 
     199        
     200      ENDIF  ! ln_bias .AND. ln_bias_pc_app 
     201 
     202 
    125203      IF( l_trddyn ) THEN      ! save the hydrostatic pressure gradient trends for momentum trend diagnostics 
    126204         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
     
    132210      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' hpg  - Ua: ', mask1=umask,   & 
    133211         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    134       ! 
    135       IF ( ln_bias .AND. ln_bias_pc_app )  THEN 
    136          IF(lwp) THEN  
    137          WRITE(numout,*) " ! restore original density" 
    138          ENDIF 
    139          rhd(:,:,:) = z_rhd_st(:,:,:)     ! restore original density 
    140          gru(:,:)   = z_gru_st(:,:) 
    141          grv(:,:)   = z_grv_st(:,:) 
    142  
    143          !Deallocate tempory variables 
    144          DEALLOCATE( z_rhd_st,     & 
    145             &        z_gru_st,     & 
    146             &        z_grv_st      ) 
    147       ENDIF 
    148212      ! 
    149213      IF( nn_timing == 1 )  CALL timing_stop('dyn_hpg') 
Note: See TracChangeset for help on using the changeset viewer.