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 11442 for branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 – NEMO

Ignore:
Timestamp:
2019-08-16T12:32:43+02:00 (5 years ago)
Author:
mattmartin
Message:

Introduction of stochastic physics in NEMO, based on Andrea Storto's code.
For details, see ticket https://code.metoffice.gov.uk/trac/utils/ticket/251.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

    r10302 r11442  
    5757   USE agrif_opa_update 
    5858#endif 
    59  
    60  
     59   USE stopack 
    6160 
    6261   IMPLICIT NONE 
     
    9493 
    9594   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau           ! depth of tke penetration (nn_htau) 
    96    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_niw          !: TKE budget- near-inertial waves term 
    97    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   efr            ! surface boundary condition for nn_etau = 4 
    9895   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl          ! now mixing lenght of dissipation 
    9996#if defined key_c1d 
     
    10299   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_pdl, e_ric   !: prandl and local Richardson numbers 
    103100#endif 
    104  
     101   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rn_lc0, rn_ediff0, rn_ediss0, rn_ebb0, rn_efr0 
     102   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_niw          !: TKE budget- near-inertial waves term 
     103   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   efr            ! surface boundary condition for nn_etau = 4 
     104    
    105105   !! * Substitutions 
    106106#  include "domzgr_substitute.h90" 
     
    184184         avmu(:,:,:) = avmu_k(:,:,:)  
    185185         avmv(:,:,:) = avmv_k(:,:,:)  
    186       ENDIF  
     186      ENDIF 
     187      ! 
     188      IF( ln_stopack) THEN 
     189         IF( nn_spp_tkelc > 0 ) THEN 
     190             rn_lc0 = rn_lc 
     191             CALL spp_gen( kt, rn_lc0, nn_spp_tkelc, rn_tkelc_sd, jk_spp_tkelc ) 
     192         ENDIF 
     193         IF( nn_spp_tkedf > 0 ) THEN 
     194             rn_ediff0 = rn_ediff 
     195             CALL spp_gen( kt, rn_ediff0, nn_spp_tkedf, rn_tkedf_sd, jk_spp_tkedf ) 
     196         ENDIF 
     197         IF( nn_spp_tkeds > 0 ) THEN 
     198             rn_ediss0 = rn_ediss 
     199             CALL spp_gen( kt, rn_ediss0, nn_spp_tkeds, rn_tkeds_sd, jk_spp_tkeds ) 
     200         ENDIF 
     201         IF( nn_spp_tkebb > 0 ) THEN 
     202             rn_ebb0 = rn_ebb 
     203             CALL spp_gen( kt, rn_ebb0, nn_spp_tkebb, rn_tkebb_sd, jk_spp_tkebb ) 
     204        ENDIF 
     205         IF( nn_spp_tkefr > 0 ) THEN 
     206             rn_efr0 = rn_efr 
     207             CALL spp_gen( kt, rn_efr0, nn_spp_tkefr, rn_tkefr_sd, jk_spp_tkefr ) 
     208         ENDIF 
     209      ENDIF 
    187210      ! 
    188211      CALL tke_tke      ! now tke (en) 
     
    199222      IF( .NOT.Agrif_Root() )   CALL Agrif_Update_Tke( kt )      ! children only 
    200223#endif       
    201      !  
     224      IF (  kt == nitend ) THEN 
     225        DEALLOCATE ( rn_lc0, rn_ediff0, rn_ediss0, rn_ebb0, rn_efr0 ) 
     226      ENDIF 
     227      ! 
     228 
    202229   END SUBROUTINE zdf_tke 
    203230 
     
    225252      REAL(wp) ::   zrhoa  = 1.22                   ! Air density kg/m3 
    226253      REAL(wp) ::   zcdrag = 1.5e-3                 ! drag coefficient 
    227       REAL(wp) ::   zbbrau, zesh2                   ! temporary scalars 
    228       REAL(wp) ::   zfact1, zfact2, zfact3          !    -         - 
     254      REAL(wp) ::   zesh2                           ! temporary scalars 
     255      REAL(wp) ::   zfact1                          !    -         - 
    229256      REAL(wp) ::   ztx2  , zty2  , zcof            !    -         - 
    230257      REAL(wp) ::   ztau  , zdif                    !    -         - 
     
    233260!!bfr      REAL(wp) ::   zebot                           !    -         - 
    234261      INTEGER , POINTER, DIMENSION(:,:  ) :: imlc 
    235       REAL(wp), POINTER, DIMENSION(:,:  ) :: zhlc 
     262      REAL(wp), POINTER, DIMENSION(:,:  ) :: zhlc, zbbrau,zfact2,zfact3 
    236263      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw 
    237264      !!-------------------------------------------------------------------- 
     
    240267      ! 
    241268      CALL wrk_alloc( jpi,jpj, imlc )    ! integer 
    242       CALL wrk_alloc( jpi,jpj, zhlc )  
     269      CALL wrk_alloc( jpi,jpj, zhlc ) 
     270      CALL wrk_alloc( jpi,jpj, zbbrau ) 
     271      CALL wrk_alloc( jpi,jpj, zfact2 ) 
     272      CALL wrk_alloc( jpi,jpj, zfact3 ) 
    243273      CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw )  
    244274      ! 
    245       zbbrau = rn_ebb / rau0       ! Local constant initialisation 
     275      zbbrau = rn_ebb0 / rau0       ! Local constant initialisation 
    246276      zfact1 = -.5_wp * rdt  
    247       zfact2 = 1.5_wp * rdt * rn_ediss 
    248       zfact3 = 0.5_wp       * rn_ediss 
     277      zfact2 = 1.5_wp * rdt * rn_ediss0 
     278      zfact3 = 0.5_wp       * rn_ediss0 
    249279      ! 
    250280      ! 
     
    261291      DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0) 
    262292         DO ji = fs_2, fs_jpim1   ! vector opt. 
    263             en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
     293            en(ji,jj,1) = MAX( rn_emin0, zbbrau(ji,jj) * taum(ji,jj) ) * tmask(ji,jj,1) 
    264294         END DO 
    265295      END DO 
     
    326356                  !                                           ! vertical velocity due to LC 
    327357                  zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) ) 
    328                   zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) ) 
     358                  zwlc = zind * rn_lc0(ji,jj) * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) ) 
    329359                  !                                           ! TKE Langmuir circulation source term 
    330                   en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( 1._wp - fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc ) /   & 
     360                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( 1._wp - fr_i(ji,jj) )* ( zwlc * zwlc * zwlc ) /   & 
    331361                     &   zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     362 
    332363               END DO 
    333364            END DO 
     
    380411               zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw) 
    381412               zd_lw(ji,jj,jk) = zzd_lw 
    382                zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk) 
     413               zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2(ji,jj) * dissl(ji,jj,jk) * tmask(ji,jj,jk) 
    383414               ! 
    384415               !                                   ! right hand side in en 
    385416               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zesh2  -   avt(ji,jj,jk) * rn2(ji,jj,jk)    & 
    386                   &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) & 
     417                  &                                 + zfact3(ji,jj) * dissl(ji,jj,jk) * en (ji,jj,jk)  ) & 
    387418                  &                                 * wmask(ji,jj,jk) 
    388419            END DO 
     
    455486            DO jj = 2, jpjm1 
    456487               DO ji = fs_2, fs_jpim1   ! vector opt. 
    457                   en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
     488                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr0(ji,jj) * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    458489                     &                                 * ( 1._wp - fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    459490               END DO 
     
    464495            DO ji = fs_2, fs_jpim1   ! vector opt. 
    465496               jk = nmln(ji,jj) 
    466                en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
     497               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr0(ji,jj) * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    467498                  &                                 * ( 1._wp - fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    468499            END DO 
     
    480511                  zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean  
    481512                  zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications... 
    482                   en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
     513                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau(ji,jj) * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    483514                     &                        * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    484515               END DO 
     
    516547      CALL wrk_dealloc( jpi,jpj, imlc )    ! integer 
    517548      CALL wrk_dealloc( jpi,jpj, zhlc )  
     549      CALL wrk_dealloc( jpi,jpj, zbbrau )  
     550      CALL wrk_dealloc( jpi,jpj, zfact2 )  
     551      CALL wrk_dealloc( jpi,jpj, zfact3 )  
    518552      CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw )  
    519553      ! 
     
    698732            DO ji = fs_2, fs_jpim1   ! vector opt. 
    699733               zsqen = SQRT( en(ji,jj,jk) ) 
    700                zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen 
     734               zav   = rn_ediff0(ji,jj) * zmxlm(ji,jj,jk) * zsqen 
    701735               avm  (ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk) 
    702736               avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 
     
    704738            END DO 
    705739         END DO 
     740         IF( ln_stopack) THEN 
     741            IF(nn_spp_avt > 0 ) CALL spp_gen( 1, avt(:,:,jk), nn_spp_avt, rn_avt_sd, jk_spp_avt, jk) 
     742            IF(nn_spp_avm > 0 ) CALL spp_gen( 1, avm(:,:,jk), nn_spp_avm, rn_avm_sd, jk_spp_avm, jk) 
     743         ENDIF 
    706744      END DO 
    707745      CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged) 
     
    840878         rn_mxl0 = rmxl_min 
    841879      ENDIF 
     880 
     881      ALLOCATE( rn_lc0   (jpi,jpj) ) ; rn_lc0    = rn_lc 
     882      ALLOCATE( rn_ediff0(jpi,jpj) ) ; rn_ediff0 = rn_ediff 
     883      ALLOCATE( rn_ediss0(jpi,jpj) ) ; rn_ediss0 = rn_ediss 
     884      ALLOCATE( rn_ebb0  (jpi,jpj) ) ; rn_ebb0   = rn_ebb 
     885      ALLOCATE( rn_efr0  (jpi,jpj) ) ; rn_efr0   = rn_efr 
    842886       
    843887      IF( nn_etau == 2  ) THEN 
     
    952996              CALL tke_avn                               ! recompute avt, avm, avmu, avmv and dissl (approximation) 
    953997              ! 
    954               avt_k (:,:,:) = avt (:,:,:) 
    955               avm_k (:,:,:) = avm (:,:,:) 
    956               avmu_k(:,:,:) = avmu(:,:,:) 
    957               avmv_k(:,:,:) = avmv(:,:,:) 
    958               ! 
    959998              DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke( jit )   ;   END DO 
    960999           ENDIF 
    9611000        ELSE                                   !* Start from rest 
    9621001           en(:,:,:) = rn_emin * tmask(:,:,:) 
    963            DO jk = 1, jpk                           ! set the Kz to the background value 
    964               avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 
    965               avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 
    966               avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk) 
    967               avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk) 
    968            END DO 
    9691002        ENDIF 
    9701003        ! 
     1004        avt_k (:,:,:) = avt (:,:,:) 
     1005        avm_k (:,:,:) = avm (:,:,:) 
     1006        avmu_k(:,:,:) = avmu(:,:,:) 
     1007        avmv_k(:,:,:) = avmv(:,:,:) 
     1008         
    9711009     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file 
    9721010        !                                   ! ------------------- 
Note: See TracChangeset for help on using the changeset viewer.