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/TRA/eosbn2.F90

    r5034 r5600  
    2222   !!             -   ! 2013-04  (F. Roquet, G. Madec)  add eos_rab, change bn2 computation and reorganize the module 
    2323   !!             -   ! 2014-09  (F. Roquet)  add TEOS-10, S-EOS, and modify EOS-80 
     24   !!             -   ! 2015-06  (P.A. Bouttier) eos_fzp functions changed to subroutines for AGRIF 
    2425   !!---------------------------------------------------------------------- 
    2526 
     
    4748   USE lbclnk         ! ocean lateral boundary conditions 
    4849   USE timing          ! Timing 
     50   USE stopar          ! Stochastic T/S fluctuations 
     51   USE stopts          ! Stochastic T/S fluctuations 
    4952 
    5053   IMPLICIT NONE 
     
    7275   PUBLIC   eos_init       ! called by istate module 
    7376 
    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 
     77   !                                !!* Namelist (nameos) * 
     78   INTEGER , PUBLIC ::   nn_eos     ! = 0/1/2 type of eq. of state and Brunt-Vaisala frequ. 
     79   LOGICAL , PUBLIC ::   ln_useCT  ! determine if eos_pt_from_ct is used to compute sst_m 
    7780 
    7881   !                                   !!!  simplified eos coefficients 
     
    313316      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pdep   ! depth                      [m] 
    314317      ! 
    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   !   -      - 
     318      INTEGER  ::   ji, jj, jk, jsmp             ! dummy loop indices 
     319      INTEGER  ::   jdof 
     320      REAL(wp) ::   zt , zh , zstemp, zs , ztm   ! local scalars 
     321      REAL(wp) ::   zn , zn0, zn1, zn2, zn3      !   -      - 
     322      REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign    ! local vectors 
    318323      !!---------------------------------------------------------------------- 
    319324      ! 
     
    324329      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    325330         ! 
    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) 
     331         ! Stochastic equation of state 
     332         IF ( ln_sto_eos ) THEN 
     333            ALLOCATE(zn0_sto(1:2*nn_sto_eos)) 
     334            ALLOCATE(zn_sto(1:2*nn_sto_eos)) 
     335            ALLOCATE(zsign(1:2*nn_sto_eos)) 
     336            DO jsmp = 1, 2*nn_sto_eos, 2 
     337              zsign(jsmp)   = 1._wp 
     338              zsign(jsmp+1) = -1._wp 
     339            END DO 
     340            ! 
     341            DO jk = 1, jpkm1 
     342               DO jj = 1, jpj 
     343                  DO ji = 1, jpi 
     344                     ! 
     345                     ! compute density (2*nn_sto_eos) times: 
     346                     ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts) 
     347                     ! (2) for t-dt, s-ds (with the opposite fluctuation) 
     348                     DO jsmp = 1, nn_sto_eos*2 
     349                        jdof   = (jsmp + 1) / 2 
     350                        zh     = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
     351                        zt     = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0    ! temperature 
     352                        zstemp = pts  (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp) 
     353                        zs     = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 )   ! square root salinity 
     354                        ztm    = tmask(ji,jj,jk)                                         ! tmask 
     355                        ! 
     356                        zn3 = EOS013*zt   & 
     357                           &   + EOS103*zs+EOS003 
     358                           ! 
     359                        zn2 = (EOS022*zt   & 
     360                           &   + EOS112*zs+EOS012)*zt   & 
     361                           &   + (EOS202*zs+EOS102)*zs+EOS002 
     362                           ! 
     363                        zn1 = (((EOS041*zt   & 
     364                           &   + EOS131*zs+EOS031)*zt   & 
     365                           &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
     366                           &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
     367                           &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
     368                           ! 
     369                        zn0_sto(jsmp) = (((((EOS060*zt   & 
     370                           &   + EOS150*zs+EOS050)*zt   & 
     371                           &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
     372                           &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
     373                           &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
     374                           &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
     375                           &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
     376                           ! 
     377                        zn_sto(jsmp)  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp) 
     378                     END DO 
     379                     ! 
     380                     ! compute stochastic density as the mean of the (2*nn_sto_eos) densities 
     381                     prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp 
     382                     DO jsmp = 1, nn_sto_eos*2 
     383                        prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp)                      ! potential density referenced at the surface 
     384                        ! 
     385                        prd(ji,jj,jk) = prd(ji,jj,jk) + (  zn_sto(jsmp) * r1_rau0 - 1._wp  )   ! density anomaly (masked) 
     386                     END DO 
     387                     prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos 
     388                     prd  (ji,jj,jk) = 0.5_wp * prd  (ji,jj,jk) * ztm / nn_sto_eos 
     389                  END DO 
    361390               END DO 
    362391            END DO 
    363          END DO 
    364          ! 
     392            DEALLOCATE(zn0_sto,zn_sto,zsign) 
     393         ! Non-stochastic equation of state 
     394         ELSE 
     395            DO jk = 1, jpkm1 
     396               DO jj = 1, jpj 
     397                  DO ji = 1, jpi 
     398                     ! 
     399                     zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
     400                     zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
     401                     zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     402                     ztm = tmask(ji,jj,jk)                                         ! tmask 
     403                     ! 
     404                     zn3 = EOS013*zt   & 
     405                        &   + EOS103*zs+EOS003 
     406                        ! 
     407                     zn2 = (EOS022*zt   & 
     408                        &   + EOS112*zs+EOS012)*zt   & 
     409                        &   + (EOS202*zs+EOS102)*zs+EOS002 
     410                        ! 
     411                     zn1 = (((EOS041*zt   & 
     412                        &   + EOS131*zs+EOS031)*zt   & 
     413                        &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
     414                        &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
     415                        &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
     416                        ! 
     417                     zn0 = (((((EOS060*zt   & 
     418                        &   + EOS150*zs+EOS050)*zt   & 
     419                        &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
     420                        &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
     421                        &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
     422                        &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
     423                        &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
     424                        ! 
     425                     zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     426                     ! 
     427                     prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface 
     428                     ! 
     429                     prd(ji,jj,jk) = (  zn * r1_rau0 - 1._wp  ) * ztm      ! density anomaly (masked) 
     430                  END DO 
     431               END DO 
     432            END DO 
     433         ENDIF 
     434          
    365435      CASE( 1 )                !==  simplified EOS  ==! 
    366436         ! 
     
    922992 
    923993 
    924    FUNCTION eos_fzp_2d( psal, pdep ) RESULT( ptf ) 
     994   SUBROUTINE  eos_fzp_2d( psal, ptf, pdep ) 
    925995      !!---------------------------------------------------------------------- 
    926996      !!                 ***  ROUTINE eos_fzp  *** 
     
    9361006      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   psal   ! salinity   [psu] 
    9371007      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ), OPTIONAL ::   pdep   ! depth      [m] 
    938       REAL(wp), DIMENSION(jpi,jpj)                          ::   ptf   ! freezing temperature [Celcius] 
     1008      REAL(wp), DIMENSION(jpi,jpj), INTENT(out  )           ::   ptf    ! freezing temperature [Celcius] 
    9391009      ! 
    9401010      INTEGER  ::   ji, jj   ! dummy loop indices 
     
    9691039         nstop = nstop + 1 
    9701040         ! 
    971       END SELECT 
    972       ! 
    973    END FUNCTION eos_fzp_2d 
    974  
    975   FUNCTION eos_fzp_0d( psal, pdep ) RESULT( ptf ) 
     1041      END SELECT       
     1042      ! 
     1043  END SUBROUTINE eos_fzp_2d 
     1044 
     1045  SUBROUTINE eos_fzp_0d( psal, ptf, pdep ) 
    9761046      !!---------------------------------------------------------------------- 
    9771047      !!                 ***  ROUTINE eos_fzp  *** 
     
    9851055      !! Reference  :   UNESCO tech. papers in the marine science no. 28. 1978 
    9861056      !!---------------------------------------------------------------------- 
    987       REAL(wp), INTENT(in)           ::   psal   ! salinity   [psu] 
    988       REAL(wp), INTENT(in), OPTIONAL ::   pdep   ! depth      [m] 
    989       REAL(wp)                       ::   ptf   ! freezing temperature [Celcius] 
     1057      REAL(wp), INTENT(in )           ::   psal         ! salinity   [psu] 
     1058      REAL(wp), INTENT(in ), OPTIONAL ::   pdep         ! depth      [m] 
     1059      REAL(wp), INTENT(out)           ::   ptf          ! freezing temperature [Celcius] 
    9901060      ! 
    9911061      REAL(wp) :: zs   ! local scalars 
     
    10171087      END SELECT 
    10181088      ! 
    1019    END FUNCTION eos_fzp_0d 
     1089   END SUBROUTINE eos_fzp_0d 
    10201090 
    10211091 
     
    11831253            WRITE(numout,*) '             model uses Conservative Temperature' 
    11841254            WRITE(numout,*) '             Important: model must be initialized with CT and SA fields' 
     1255         ELSE 
     1256            WRITE(numout,*) '             model does not use Conservative Temperature' 
    11851257         ENDIF 
    11861258      ENDIF 
     
    15891661      END SELECT 
    15901662      ! 
     1663      rau0_rcp    = rau0 * rcp  
    15911664      r1_rau0     = 1._wp / rau0 
    15921665      r1_rcp      = 1._wp / rcp 
    1593       r1_rau0_rcp = 1._wp / ( rau0 * rcp ) 
     1666      r1_rau0_rcp = 1._wp / rau0_rcp  
    15941667      ! 
    15951668      IF(lwp) WRITE(numout,*) 
     
    15971670      IF(lwp) WRITE(numout,*) '          1. / rau0                        r1_rau0  = ', r1_rau0, ' m^3/kg' 
    15981671      IF(lwp) WRITE(numout,*) '          ocean specific heat                 rcp   = ', rcp    , ' J/Kelvin' 
     1672      IF(lwp) WRITE(numout,*) '          rau0 * rcp                       rau0_rcp = ', rau0_rcp 
    15991673      IF(lwp) WRITE(numout,*) '          1. / ( rau0 * rcp )           r1_rau0_rcp = ', r1_rau0_rcp 
    16001674      ! 
Note: See TracChangeset for help on using the changeset viewer.