Ignore:
Timestamp:
2015-07-15T17:46:12+02:00 (5 years ago)
Author:
andrewryan
Message:

merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r5034 r5600  
    1616   !!            3.4  !  2011-11  (H. Liu) hpg_prj: Original code for s-coordinates 
    1717   !!                 !           (A. Coward) suppression of hel, wdj and rot options 
     18   !!            3.6  !  2014-11  (P. Mathiot) hpg_isf: original code for ice shelf cavity 
    1819   !!---------------------------------------------------------------------- 
    1920 
     
    2526   !!       hpg_zps  : z-coordinate plus partial steps (interpolation) 
    2627   !!       hpg_sco  : s-coordinate (standard jacobian formulation) 
     28   !!       hpg_isf  : s-coordinate (sco formulation) adapted to ice shelf 
    2729   !!       hpg_djc  : s-coordinate (Density Jacobian with Cubic polynomial) 
    2830   !!       hpg_prj  : s-coordinate (Pressure Jacobian with Cubic polynomial) 
     
    5557   LOGICAL , PUBLIC ::   ln_hpg_djc      !: s-coordinate (Density Jacobian with Cubic polynomial) 
    5658   LOGICAL , PUBLIC ::   ln_hpg_prj      !: s-coordinate (Pressure Jacobian scheme) 
     59   LOGICAL , PUBLIC ::   ln_hpg_isf      !: s-coordinate similar to sco modify for isf 
    5760   LOGICAL , PUBLIC ::   ln_dynhpg_imp   !: semi-implicite hpg flag 
    5861 
     
    97100      CASE (  3 )   ;   CALL hpg_djc    ( kt )      ! s-coordinate (Density Jacobian with Cubic polynomial) 
    98101      CASE (  4 )   ;   CALL hpg_prj    ( kt )      ! s-coordinate (Pressure Jacobian scheme) 
     102      CASE (  5 )   ;   CALL hpg_isf    ( kt )      ! s-coordinate similar to sco modify for ice shelf 
    99103      END SELECT 
    100104      ! 
     
    128132      !! 
    129133      NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco,     & 
    130          &                 ln_hpg_djc, ln_hpg_prj, ln_dynhpg_imp 
     134         &                 ln_hpg_djc, ln_hpg_prj, ln_hpg_isf, ln_dynhpg_imp 
    131135      !!---------------------------------------------------------------------- 
    132136      ! 
     
    148152         WRITE(numout,*) '      z-coord. - partial steps (interpolation)          ln_hpg_zps    = ', ln_hpg_zps 
    149153         WRITE(numout,*) '      s-coord. (standard jacobian formulation)          ln_hpg_sco    = ', ln_hpg_sco 
     154         WRITE(numout,*) '      s-coord. (standard jacobian formulation) for isf  ln_hpg_isf    = ', ln_hpg_isf 
    150155         WRITE(numout,*) '      s-coord. (Density Jacobian: Cubic polynomial)     ln_hpg_djc    = ', ln_hpg_djc 
    151156         WRITE(numout,*) '      s-coord. (Pressure Jacobian: Cubic polynomial)    ln_hpg_prj    = ', ln_hpg_prj 
     
    158163                           & either  ln_hpg_sco or  ln_hpg_prj instead') 
    159164      ! 
    160       IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj) )   & 
     165      IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) )   & 
    161166         &   CALL ctl_stop('dyn_hpg_init : variable volume key_vvl requires:& 
    162167                           & the standard jacobian formulation hpg_sco or & 
    163168                           & the pressure jacobian formulation hpg_prj') 
     169 
     170      IF(       ln_hpg_isf .AND. .NOT. ln_isfcav )   & 
     171         &   CALL ctl_stop( ' hpg_isf not available if ln_isfcav = false ' ) 
     172      IF( .NOT. ln_hpg_isf .AND.       ln_isfcav )   & 
     173         &   CALL ctl_stop( 'Only hpg_isf has been corrected to work with ice shelf cavity.' ) 
    164174      ! 
    165175      !                               ! Set nhpg from ln_hpg_... flags 
     
    169179      IF( ln_hpg_djc )   nhpg = 3 
    170180      IF( ln_hpg_prj )   nhpg = 4 
     181      IF( ln_hpg_isf )   nhpg = 5 
    171182      ! 
    172183      !                               ! Consistency check 
     
    177188      IF( ln_hpg_djc )   ioptio = ioptio + 1 
    178189      IF( ln_hpg_prj )   ioptio = ioptio + 1 
     190      IF( ln_hpg_isf )   ioptio = ioptio + 1 
    179191      IF( ioptio /= 1 )   CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 
    180       IF( (ln_hpg_zco .OR. ln_hpg_zps .OR. ln_hpg_djc .OR. ln_hpg_prj ) .AND. nn_isf .NE. 0 )   & 
    181           &  CALL ctl_stop( 'Only hpg_sco has been corrected to work with ice shelf cavity.' ) 
     192      !  
     193      ! initialisation of ice load 
     194      riceload(:,:)=0.0 
    182195      ! 
    183196   END SUBROUTINE dyn_hpg_init 
     
    345358   END SUBROUTINE hpg_zps 
    346359 
    347  
    348360   SUBROUTINE hpg_sco( kt ) 
    349361      !!--------------------------------------------------------------------- 
     
    366378      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    367379      !! 
     380      INTEGER  ::   ji, jj, jk                 ! dummy loop indices 
     381      REAL(wp) ::   zcoef0, zuap, zvap, znad   ! temporary scalars 
     382      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj 
     383      !!---------------------------------------------------------------------- 
     384      ! 
     385      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 
     386      ! 
     387      IF( kt == nit000 ) THEN 
     388         IF(lwp) WRITE(numout,*) 
     389         IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend' 
     390         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, OPA original scheme used' 
     391      ENDIF 
     392 
     393      ! Local constant initialization 
     394      zcoef0 = - grav * 0.5_wp 
     395      ! To use density and not density anomaly 
     396      IF ( lk_vvl ) THEN   ;     znad = 1._wp          ! Variable volume 
     397      ELSE                 ;     znad = 0._wp         ! Fixed volume 
     398      ENDIF 
     399 
     400      ! Surface value 
     401      DO jj = 2, jpjm1 
     402         DO ji = fs_2, fs_jpim1   ! vector opt. 
     403            ! hydrostatic pressure gradient along s-surfaces 
     404            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj  ,1) * ( znad + rhd(ji+1,jj  ,1) )   & 
     405               &                                  - fse3w(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) ) 
     406            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji  ,jj+1,1) * ( znad + rhd(ji  ,jj+1,1) )   & 
     407               &                                  - fse3w(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) ) 
     408            ! s-coordinate pressure gradient correction 
     409            zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) + 2._wp * znad )   & 
     410               &           * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) 
     411            zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) + 2._wp * znad )   & 
     412               &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 
     413            ! add to the general momentum trend 
     414            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 
     415            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap 
     416         END DO 
     417      END DO 
     418 
     419      ! interior value (2=<jk=<jpkm1) 
     420      DO jk = 2, jpkm1 
     421         DO jj = 2, jpjm1 
     422            DO ji = fs_2, fs_jpim1   ! vector opt. 
     423               ! hydrostatic pressure gradient along s-surfaces 
     424               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   & 
     425                  &           * (  fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   & 
     426                  &              - fse3w(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  ) 
     427               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   & 
     428                  &           * (  fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   & 
     429                  &              - fse3w(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  ) 
     430               ! s-coordinate pressure gradient correction 
     431               zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
     432                  &           * ( fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) 
     433               zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
     434                  &           * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 
     435               ! add to the general momentum trend 
     436               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 
     437               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap 
     438            END DO 
     439         END DO 
     440      END DO 
     441      ! 
     442      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 
     443      ! 
     444   END SUBROUTINE hpg_sco 
     445 
     446   SUBROUTINE hpg_isf( kt ) 
     447      !!--------------------------------------------------------------------- 
     448      !!                  ***  ROUTINE hpg_sco  *** 
     449      !! 
     450      !! ** Method  :   s-coordinate case. Jacobian scheme. 
     451      !!      The now hydrostatic pressure gradient at a given level, jk, 
     452      !!      is computed by taking the vertical integral of the in-situ 
     453      !!      density gradient along the model level from the suface to that 
     454      !!      level. s-coordinates (ln_sco): a corrective term is added 
     455      !!      to the horizontal pressure gradient : 
     456      !!         zhpi = grav .....  + 1/e1u mi(rhd) di[ grav dep3w ] 
     457      !!         zhpj = grav .....  + 1/e2v mj(rhd) dj[ grav dep3w ] 
     458      !!      add it to the general momentum trend (ua,va). 
     459      !!         ua = ua - 1/e1u * zhpi 
     460      !!         va = va - 1/e2v * zhpj 
     461      !!      iceload is added and partial cell case are added to the top and bottom 
     462      !!       
     463      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 
     464      !!---------------------------------------------------------------------- 
     465      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
     466      !! 
    368467      INTEGER  ::   ji, jj, jk, iku, ikv, ikt, iktp1i, iktp1j                 ! dummy loop indices 
    369468      REAL(wp) ::   zcoef0, zuap, zvap, znad, ze3wu, ze3wv, zuapint, zvapint, zhpjint, zhpiint, zdzwt, zdzwtjp1, zdzwtip1             ! temporary scalars 
     
    379478     IF( kt == nit000 ) THEN 
    380479         IF(lwp) WRITE(numout,*) 
    381          IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend' 
     480         IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf' 
    382481         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, OPA original scheme used' 
    383482      ENDIF 
     
    565664!================================================================================== 
    566665 
    567 # if defined key_vectopt_loop 
    568          jj = 1 
    569          DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
    570 # else 
    571666      DO jj = 2, jpjm1 
    572667         DO ji = 2, jpim1 
    573 # endif 
    574668            iku = mbku(ji,jj) 
    575669            ikv = mbkv(ji,jj) 
     
    598692               va(ji,jj,ikv)   =  va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap 
    599693            END IF 
    600 # if ! defined key_vectopt_loop 
    601          END DO 
    602 # endif 
     694         END DO 
    603695      END DO 
    604696      
     
    610702      CALL wrk_dealloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 
    611703      ! 
    612    END SUBROUTINE hpg_sco 
     704   END SUBROUTINE hpg_isf 
    613705 
    614706 
     
    864956      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdept, zrhh 
    865957      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 
     958      REAL(wp), POINTER, DIMENSION(:,:)   ::   zsshu_n, zsshv_n 
    866959      !!---------------------------------------------------------------------- 
    867960      ! 
    868961      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 
    869962      CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 
     963      CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n ) 
    870964      ! 
    871965      IF( kt == nit000 ) THEN 
     
    9481042 
    9491043      ! Z coordinate of U(ji,jj,1:jpkm1) and V(ji,jj,1:jpkm1) 
     1044 
     1045      ! Prepare zsshu_n and zsshv_n 
    9501046      DO jj = 2, jpjm1 
    9511047        DO ji = 2, jpim1 
    952           zu(ji,jj,1) = - ( fse3u(ji,jj,1) - sshn(ji,jj) * znad)    ! probable bug: changed from sshu_n for ztilde compilation 
    953           zv(ji,jj,1) = - ( fse3v(ji,jj,1) - sshn(ji,jj) * znad)    ! probable bug: changed from sshv_n for ztilde compilation 
     1048          zsshu_n(ji,jj) = (e12u(ji,jj) * sshn(ji,jj) + e12u(ji+1, jj) * sshn(ji+1,jj)) * & 
     1049                         & r1_e12u(ji,jj) * umask(ji,jj,1) * 0.5_wp  
     1050          zsshv_n(ji,jj) = (e12v(ji,jj) * sshn(ji,jj) + e12v(ji+1, jj) * sshn(ji,jj+1)) * & 
     1051                         & r1_e12v(ji,jj) * vmask(ji,jj,1) * 0.5_wp  
     1052        END DO 
     1053      END DO 
     1054 
     1055      DO jj = 2, jpjm1 
     1056        DO ji = 2, jpim1 
     1057          zu(ji,jj,1) = - ( fse3u(ji,jj,1) - zsshu_n(ji,jj) * znad)  
     1058          zv(ji,jj,1) = - ( fse3v(ji,jj,1) - zsshv_n(ji,jj) * znad) 
    9541059        END DO 
    9551060      END DO 
     
    11131218      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 
    11141219      CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 
     1220      CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n ) 
    11151221      ! 
    11161222   END SUBROUTINE hpg_prj 
Note: See TracChangeset for help on using the changeset viewer.