Ignore:
Timestamp:
2015-05-27T12:47:55+02:00 (6 years ago)
Author:
pabouttier
Message:

Commit stochastic parametrization module and perturbation of EOS

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5177_CNRS4_stopar/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90

    r5147 r5296  
    3030   !!   eos_insitu_2d  : Compute the in situ density for 2d fields 
    3131   !!   bn2            : Compute the Brunt-Vaisala frequency 
    32    !!   eos_rab        : generic interface of in situ thermal/haline expansion ratio  
     32   !!   eos_rab        : generic interface of in situ thermal/haline expansion ratio 
    3333   !!   eos_rab_3d     : compute in situ thermal/haline expansion ratio 
    3434   !!   eos_rab_2d     : compute in situ thermal/haline expansion ratio for 2d fields 
     
    4242   USE in_out_manager  ! I/O manager 
    4343   USE lib_mpp         ! MPP library 
    44    USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     44   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
    4545   USE prtctl          ! Print control 
    4646   USE wrk_nemo        ! Memory Allocation 
    47    USE lbclnk         ! ocean lateral boundary conditions 
     47   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 
     
    6062   END INTERFACE 
    6163   ! 
    62    INTERFACE eos_fzp  
     64   INTERFACE eos_fzp 
    6365      MODULE PROCEDURE eos_fzp_2d, eos_fzp_0d 
    6466   END INTERFACE 
     
    7880   !                                   !!!  simplified eos coefficients 
    7981   ! default value: Vallis 2006 
    80    REAL(wp) ::   rn_a0      = 1.6550e-1_wp     ! thermal expansion coeff.  
    81    REAL(wp) ::   rn_b0      = 7.6554e-1_wp     ! saline  expansion coeff.  
    82    REAL(wp) ::   rn_lambda1 = 5.9520e-2_wp     ! cabbeling coeff. in T^2         
    83    REAL(wp) ::   rn_lambda2 = 5.4914e-4_wp     ! cabbeling coeff. in S^2         
    84    REAL(wp) ::   rn_mu1     = 1.4970e-4_wp     ! thermobaric coeff. in T   
    85    REAL(wp) ::   rn_mu2     = 1.1090e-5_wp     ! thermobaric coeff. in S   
    86    REAL(wp) ::   rn_nu      = 2.4341e-3_wp     ! cabbeling coeff. in theta*salt   
    87     
     82   REAL(wp) ::   rn_a0      = 1.6550e-1_wp     ! thermal expansion coeff. 
     83   REAL(wp) ::   rn_b0      = 7.6554e-1_wp     ! saline  expansion coeff. 
     84   REAL(wp) ::   rn_lambda1 = 5.9520e-2_wp     ! cabbeling coeff. in T^2 
     85   REAL(wp) ::   rn_lambda2 = 5.4914e-4_wp     ! cabbeling coeff. in S^2 
     86   REAL(wp) ::   rn_mu1     = 1.4970e-4_wp     ! thermobaric coeff. in T 
     87   REAL(wp) ::   rn_mu2     = 1.1090e-5_wp     ! thermobaric coeff. in S 
     88   REAL(wp) ::   rn_nu      = 2.4341e-3_wp     ! cabbeling coeff. in theta*salt 
     89 
    8890   ! TEOS10/EOS80 parameters 
    8991   REAL(wp) ::   r1_S0, r1_T0, r1_Z0, rdeltaS 
    90     
     92 
    9193   ! EOS parameters 
    9294   REAL(wp) ::   EOS000 , EOS100 , EOS200 , EOS300 , EOS400 , EOS500 , EOS600 
     
    106108   REAL(wp) ::   EOS022 
    107109   REAL(wp) ::   EOS003 , EOS103 
    108    REAL(wp) ::   EOS013  
    109     
     110   REAL(wp) ::   EOS013 
     111 
    110112   ! ALPHA parameters 
    111113   REAL(wp) ::   ALP000 , ALP100 , ALP200 , ALP300 , ALP400 , ALP500 
     
    122124   REAL(wp) ::   ALP012 
    123125   REAL(wp) ::   ALP003 
    124     
     126 
    125127   ! BETA parameters 
    126128   REAL(wp) ::   BET000 , BET100 , BET200 , BET300 , BET400 , BET500 
     
    149151   REAL(wp) ::   PEN002 , PEN102 
    150152   REAL(wp) ::   PEN012 
    151     
     153 
    152154   ! ALPHA_PEN parameters 
    153155   REAL(wp) ::   APE000 , APE100 , APE200 , APE300 
     
    279281                     &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   & 
    280282                     &  - rn_nu * zt * zs 
    281                      !                                  
     283                     ! 
    282284                  prd(ji,jj,jk) = zn * r1_rau0 * ztm                ! density anomaly (masked) 
    283285               END DO 
     
    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, jdof       ! dummy loop indices 
     318      REAL(wp) ::   zt , zh , zstemp, zs , ztm   ! local scalars 
     319      REAL(wp) ::   zn , zn0, zn1, zn2, zn3      !   -      - 
     320      REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign            ! local vectors 
    318321      !!---------------------------------------------------------------------- 
    319322      ! 
     
    324327      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    325328         ! 
    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) 
     329         ! Stochastic equation of state 
     330         IF ( ln_sto_eos ) THEN 
     331            ALLOCATE(zn0_sto(1:2*nn_sto_eos)) 
     332            ALLOCATE(zn_sto(1:2*nn_sto_eos)) 
     333            ALLOCATE(zsign(1:2*nn_sto_eos)) 
     334            DO jsmp = 1, 2*nn_sto_eos, 2 
     335              zsign(jsmp)   = 1._wp 
     336              zsign(jsmp+1) = -1._wp 
     337            END DO 
     338            ! 
     339            DO jk = 1, jpkm1 
     340               DO jj = 1, jpj 
     341                  DO ji = 1, jpi 
     342                     ! 
     343                     ! compute density (2*nn_sto_eos) times: 
     344                     ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts) 
     345                     ! (2) for t-dt, s-ds (with the opposite fluctuation) 
     346                     DO jsmp = 1, nn_sto_eos*2 
     347                        zh     = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
     348                        zt     = pts (ji,jj,jk,jp_tem) * r1_T0  + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)  ! temperature 
     349                        zstemp = pts   (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp) 
     350                        zs     = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 )   ! square root salinity 
     351                        ztm    = tmask(ji,jj,jk)                                         ! tmask 
     352                        ! 
     353                        zn3 = EOS013*zt   & 
     354                           &   + EOS103*zs+EOS003 
     355                           ! 
     356                        zn2 = (EOS022*zt   & 
     357                           &   + EOS112*zs+EOS012)*zt   & 
     358                           &   + (EOS202*zs+EOS102)*zs+EOS002 
     359                           ! 
     360                        zn1 = (((EOS041*zt   & 
     361                           &   + EOS131*zs+EOS031)*zt   & 
     362                           &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
     363                           &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
     364                           &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
     365                           ! 
     366                        zn0_sto(jsmp) = (((((EOS060*zt   & 
     367                           &   + EOS150*zs+EOS050)*zt   & 
     368                           &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
     369                           &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
     370                           &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
     371                           &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
     372                           &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
     373                           ! 
     374                        zn_sto(jsmp)  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     375                     END DO 
     376                     ! 
     377                     ! 
     378                     ! compute stochastic density as the mean of the (2*nn_sto_eos) densities 
     379                     prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp 
     380                     DO jsmp = 1, nn_sto_eos*2 
     381                        prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp)                      ! potential density referenced at the surface 
     382                        !  
     383                        prd(ji,jj,jk) = prd(ji,jj,jk) + (  zn_sto(jsmp) * r1_rau0 - 1._wp  )   ! density anomaly (masked) 
     384                     END DO 
     385                     prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos 
     386                     prd  (ji,jj,jk) = 0.5_wp * prd  (ji,jj,jk) * ztm / nn_sto_eos 
     387                  END DO 
    361388               END DO 
    362389            END DO 
    363          END DO 
     390            DEALLOCATE(zn0_sto,zn_sto,zsign) 
     391         ! Non-stochastic equation of state 
     392         ELSE 
     393            DO jk = 1, jpkm1 
     394               DO jj = 1, jpj 
     395                  DO ji = 1, jpi 
     396                     ! 
     397                     zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
     398                     zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
     399                     zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     400                     ztm = tmask(ji,jj,jk)                                         ! tmask 
     401                     ! 
     402                     zn3 = EOS013*zt   & 
     403                        &   + EOS103*zs+EOS003 
     404                        ! 
     405                     zn2 = (EOS022*zt   & 
     406                        &   + EOS112*zs+EOS012)*zt   & 
     407                        &   + (EOS202*zs+EOS102)*zs+EOS002 
     408                        ! 
     409                     zn1 = (((EOS041*zt   & 
     410                        &   + EOS131*zs+EOS031)*zt   & 
     411                        &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
     412                        &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
     413                        &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
     414                        ! 
     415                     zn0 = (((((EOS060*zt   & 
     416                        &   + EOS150*zs+EOS050)*zt   & 
     417                        &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
     418                        &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
     419                        &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
     420                        &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
     421                        &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
     422                        ! 
     423                     zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     424                     ! 
     425                     prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface 
     426                     ! 
     427                     prd(ji,jj,jk) = (  zn * r1_rau0 - 1._wp  ) * ztm      ! density anomaly (masked) 
     428                  END DO 
     429               END DO 
     430            END DO 
     431         ENDIF 
    364432         ! 
    365433      CASE( 1 )                !==  simplified EOS  ==! 
     
    681749         ! 
    682750         CALL lbc_lnk( pab(:,:,jp_tem), 'T', 1. )                    ! Lateral boundary conditions 
    683          CALL lbc_lnk( pab(:,:,jp_sal), 'T', 1. )                     
     751         CALL lbc_lnk( pab(:,:,jp_sal), 'T', 1. ) 
    684752         ! 
    685753      CASE( 1 )                  !==  simplified EOS  ==! 
     
    702770         ! 
    703771         CALL lbc_lnk( pab(:,:,jp_tem), 'T', 1. )                    ! Lateral boundary conditions 
    704          CALL lbc_lnk( pab(:,:,jp_sal), 'T', 1. )                     
     772         CALL lbc_lnk( pab(:,:,jp_sal), 'T', 1. ) 
    705773         ! 
    706774      CASE DEFAULT 
     
    820888      !!                  ***  ROUTINE bn2  *** 
    821889      !! 
    822       !! ** Purpose :   Compute the local Brunt-Vaisala frequency at the  
     890      !! ** Purpose :   Compute the local Brunt-Vaisala frequency at the 
    823891      !!                time-step of the input arguments 
    824892      !! 
     
    827895      !!      N.B. N^2 is set one for all to zero at jk=1 in istate module. 
    828896      !! 
    829       !! ** Action  :   pn2 : square of the brunt-vaisala frequency at w-point  
     897      !! ** Action  :   pn2 : square of the brunt-vaisala frequency at w-point 
    830898      !! 
    831899      !!---------------------------------------------------------------------- 
     
    844912            DO ji = 1, jpi 
    845913               zrw =   ( fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk) )   & 
    846                   &  / ( fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk) )  
    847                   ! 
    848                zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw  
     914                  &  / ( fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk) ) 
     915                  ! 
     916               zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw 
    849917               zbw = pab(ji,jj,jk,jp_sal) * (1. - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw 
    850918               ! 
     
    10261094      !! ** Purpose :   Calculates nonlinear anomalies of alpha_PE, beta_PE and PE at T-points 
    10271095      !! 
    1028       !! ** Method  :   PE is defined analytically as the vertical  
     1096      !! ** Method  :   PE is defined analytically as the vertical 
    10291097      !!                   primitive of EOS times -g integrated between 0 and z>0. 
    10301098      !!                pen is the nonlinear bsq-PE anomaly: pen = ( PE - rau0 gz ) / rau0 gz - rd 
    1031       !!                                                      = 1/z * /int_0^z rd dz - rd  
     1099      !!                                                      = 1/z * /int_0^z rd dz - rd 
    10321100      !!                                where rd is the density anomaly (see eos_rhd function) 
    10331101      !!                ab_pe are partial derivatives of PE anomaly with respect to T and S: 
     
    10941162                     ! 
    10951163                  zn  = ( zn2 * zh + zn1 ) * zh + zn0 
    1096                   !                               
     1164                  ! 
    10971165                  pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rau0 * ztm 
    10981166                  ! 
     
    11091177                     ! 
    11101178                  zn  = ( zn2 * zh + zn1 ) * zh + zn0 
    1111                   !                               
     1179                  ! 
    11121180                  pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rau0 * ztm 
    11131181                  ! 
     
    15891657      END SELECT 
    15901658      ! 
    1591       rau0_rcp    = rau0 * rcp  
     1659      rau0_rcp    = rau0 * rcp 
    15921660      r1_rau0     = 1._wp / rau0 
    15931661      r1_rcp      = 1._wp / rcp 
    1594       r1_rau0_rcp = 1._wp / rau0_rcp  
     1662      r1_rau0_rcp = 1._wp / rau0_rcp 
    15951663      ! 
    15961664      IF(lwp) WRITE(numout,*) 
Note: See TracChangeset for help on using the changeset viewer.