New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 5447 for branches/UKMO/dev_r5107_mld_zint/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90 – NEMO

Ignore:
Timestamp:
2015-06-19T18:07:11+02:00 (9 years ago)
Author:
davestorkey
Message:

Update UKMO/dev_r5107_mld_zint branch to revision 5442 of the trunk.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5107_mld_zint/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90

    r5247 r5447  
    4747   USE lbclnk         ! ocean lateral boundary conditions 
    4848   USE timing          ! Timing 
     49   USE stopar          ! Stochastic T/S fluctuations 
     50   USE stopts          ! Stochastic T/S fluctuations 
    4951 
    5052   IMPLICIT NONE 
     
    7274   PUBLIC   eos_init       ! called by istate module 
    7375 
    74    !                                          !!* Namelist (nameos) * 
    75    INTEGER , PUBLIC ::   nn_eos   = 0         !: = 0/1/2 type of eq. of state and Brunt-Vaisala frequ. 
    76    LOGICAL , PUBLIC ::   ln_useCT  = .FALSE. ! determine if eos_pt_from_ct is used to compute sst_m 
     76   !                                !!* Namelist (nameos) * 
     77   INTEGER , PUBLIC ::   nn_eos     ! = 0/1/2 type of eq. of state and Brunt-Vaisala frequ. 
     78   LOGICAL , PUBLIC ::   ln_useCT  ! determine if eos_pt_from_ct is used to compute sst_m 
    7779 
    7880   !                                   !!!  simplified eos coefficients 
     
    313315      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pdep   ! depth                      [m] 
    314316      ! 
    315       INTEGER  ::   ji, jj, jk                ! dummy loop indices 
    316       REAL(wp) ::   zt , zh , zs , ztm        ! local scalars 
    317       REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      - 
     317      INTEGER  ::   ji, jj, jk, jsmp             ! dummy loop indices 
     318      INTEGER  ::   jdof 
     319      REAL(wp) ::   zt , zh , zstemp, zs , ztm   ! local scalars 
     320      REAL(wp) ::   zn , zn0, zn1, zn2, zn3      !   -      - 
     321      REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign    ! local vectors 
    318322      !!---------------------------------------------------------------------- 
    319323      ! 
     
    324328      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    325329         ! 
    326          DO jk = 1, jpkm1 
    327             DO jj = 1, jpj 
    328                DO ji = 1, jpi 
    329                   ! 
    330                   zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
    331                   zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
    332                   zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
    333                   ztm = tmask(ji,jj,jk)                                         ! tmask 
    334                   ! 
    335                   zn3 = EOS013*zt   & 
    336                      &   + EOS103*zs+EOS003 
    337                      ! 
    338                   zn2 = (EOS022*zt   & 
    339                      &   + EOS112*zs+EOS012)*zt   & 
    340                      &   + (EOS202*zs+EOS102)*zs+EOS002 
    341                      ! 
    342                   zn1 = (((EOS041*zt   & 
    343                      &   + EOS131*zs+EOS031)*zt   & 
    344                      &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
    345                      &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
    346                      &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
    347                      ! 
    348                   zn0 = (((((EOS060*zt   & 
    349                      &   + EOS150*zs+EOS050)*zt   & 
    350                      &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
    351                      &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
    352                      &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
    353                      &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
    354                      &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
    355                      ! 
    356                   zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    357                   ! 
    358                   prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface 
    359                   ! 
    360                   prd(ji,jj,jk) = (  zn * r1_rau0 - 1._wp  ) * ztm      ! density anomaly (masked) 
     330         ! Stochastic equation of state 
     331         IF ( ln_sto_eos ) THEN 
     332            ALLOCATE(zn0_sto(1:2*nn_sto_eos)) 
     333            ALLOCATE(zn_sto(1:2*nn_sto_eos)) 
     334            ALLOCATE(zsign(1:2*nn_sto_eos)) 
     335            DO jsmp = 1, 2*nn_sto_eos, 2 
     336              zsign(jsmp)   = 1._wp 
     337              zsign(jsmp+1) = -1._wp 
     338            END DO 
     339            ! 
     340            DO jk = 1, jpkm1 
     341               DO jj = 1, jpj 
     342                  DO ji = 1, jpi 
     343                     ! 
     344                     ! compute density (2*nn_sto_eos) times: 
     345                     ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts) 
     346                     ! (2) for t-dt, s-ds (with the opposite fluctuation) 
     347                     DO jsmp = 1, nn_sto_eos*2 
     348                        jdof   = (jsmp + 1) / 2 
     349                        zh     = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
     350                        zt     = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0    ! temperature 
     351                        zstemp = pts  (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp) 
     352                        zs     = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 )   ! square root salinity 
     353                        ztm    = tmask(ji,jj,jk)                                         ! tmask 
     354                        ! 
     355                        zn3 = EOS013*zt   & 
     356                           &   + EOS103*zs+EOS003 
     357                           ! 
     358                        zn2 = (EOS022*zt   & 
     359                           &   + EOS112*zs+EOS012)*zt   & 
     360                           &   + (EOS202*zs+EOS102)*zs+EOS002 
     361                           ! 
     362                        zn1 = (((EOS041*zt   & 
     363                           &   + EOS131*zs+EOS031)*zt   & 
     364                           &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
     365                           &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
     366                           &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
     367                           ! 
     368                        zn0_sto(jsmp) = (((((EOS060*zt   & 
     369                           &   + EOS150*zs+EOS050)*zt   & 
     370                           &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
     371                           &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
     372                           &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
     373                           &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
     374                           &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
     375                           ! 
     376                        zn_sto(jsmp)  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp) 
     377                     END DO 
     378                     ! 
     379                     ! compute stochastic density as the mean of the (2*nn_sto_eos) densities 
     380                     prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp 
     381                     DO jsmp = 1, nn_sto_eos*2 
     382                        prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp)                      ! potential density referenced at the surface 
     383                        ! 
     384                        prd(ji,jj,jk) = prd(ji,jj,jk) + (  zn_sto(jsmp) * r1_rau0 - 1._wp  )   ! density anomaly (masked) 
     385                     END DO 
     386                     prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos 
     387                     prd  (ji,jj,jk) = 0.5_wp * prd  (ji,jj,jk) * ztm / nn_sto_eos 
     388                  END DO 
    361389               END DO 
    362390            END DO 
    363          END DO 
    364          ! 
     391            DEALLOCATE(zn0_sto,zn_sto,zsign) 
     392         ! Non-stochastic equation of state 
     393         ELSE 
     394            DO jk = 1, jpkm1 
     395               DO jj = 1, jpj 
     396                  DO ji = 1, jpi 
     397                     ! 
     398                     zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
     399                     zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
     400                     zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     401                     ztm = tmask(ji,jj,jk)                                         ! tmask 
     402                     ! 
     403                     zn3 = EOS013*zt   & 
     404                        &   + EOS103*zs+EOS003 
     405                        ! 
     406                     zn2 = (EOS022*zt   & 
     407                        &   + EOS112*zs+EOS012)*zt   & 
     408                        &   + (EOS202*zs+EOS102)*zs+EOS002 
     409                        ! 
     410                     zn1 = (((EOS041*zt   & 
     411                        &   + EOS131*zs+EOS031)*zt   & 
     412                        &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
     413                        &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
     414                        &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
     415                        ! 
     416                     zn0 = (((((EOS060*zt   & 
     417                        &   + EOS150*zs+EOS050)*zt   & 
     418                        &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
     419                        &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
     420                        &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
     421                        &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
     422                        &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
     423                        ! 
     424                     zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     425                     ! 
     426                     prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface 
     427                     ! 
     428                     prd(ji,jj,jk) = (  zn * r1_rau0 - 1._wp  ) * ztm      ! density anomaly (masked) 
     429                  END DO 
     430               END DO 
     431            END DO 
     432         ENDIF 
     433          
    365434      CASE( 1 )                !==  simplified EOS  ==! 
    366435         ! 
     
    11831252            WRITE(numout,*) '             model uses Conservative Temperature' 
    11841253            WRITE(numout,*) '             Important: model must be initialized with CT and SA fields' 
     1254         ELSE 
     1255            WRITE(numout,*) '             model does not use Conservative Temperature' 
    11851256         ENDIF 
    11861257      ENDIF 
     
    15891660      END SELECT 
    15901661      ! 
     1662      rau0_rcp    = rau0 * rcp  
    15911663      r1_rau0     = 1._wp / rau0 
    15921664      r1_rcp      = 1._wp / rcp 
    1593       r1_rau0_rcp = 1._wp / ( rau0 * rcp ) 
     1665      r1_rau0_rcp = 1._wp / rau0_rcp  
    15941666      ! 
    15951667      IF(lwp) WRITE(numout,*) 
     
    15971669      IF(lwp) WRITE(numout,*) '          1. / rau0                        r1_rau0  = ', r1_rau0, ' m^3/kg' 
    15981670      IF(lwp) WRITE(numout,*) '          ocean specific heat                 rcp   = ', rcp    , ' J/Kelvin' 
     1671      IF(lwp) WRITE(numout,*) '          rau0 * rcp                       rau0_rcp = ', rau0_rcp 
    15991672      IF(lwp) WRITE(numout,*) '          1. / ( rau0 * rcp )           r1_rau0_rcp = ', r1_rau0_rcp 
    16001673      ! 
Note: See TracChangeset for help on using the changeset viewer.