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

Ignore:
Timestamp:
2019-07-31T18:05:50+02:00 (5 years ago)
Author:
mattmartin
Message:

Included Andrea Storto's STOPACK code into NEMO3.6 branch.

File:
1 edited

Legend:

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

    r5407 r11384  
    5353   USE timing         ! Timing 
    5454   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     55   USE stopack 
    5556 
    5657   IMPLICIT NONE 
     
    8586   REAL(wp) ::   rhftau_scl = 1.0_wp       ! scale factor applied to HF part of taum  (nn_etau=3) 
    8687 
    87    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   en             !: now turbulent kinetic energy   [m2/s2] 
     88   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rn_lc0, rn_ediff0, rn_ediss0, rn_ebb0, rn_efr0 
     89 
    8890   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau           ! depth of tke penetration (nn_htau) 
    8991   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl          ! now mixing lenght of dissipation 
    90    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avt_k , avm_k  ! not enhanced Kz 
    91    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avmu_k, avmv_k ! not enhanced Kz 
     92 
    9293#if defined key_c1d 
    9394   !                                                                        !!** 1D cfg only  **   ('key_c1d') 
     
    115116         &      e_pdl(jpi,jpj,jpk) , e_ric(jpi,jpj,jpk) ,                          & 
    116117#endif 
    117          &      en    (jpi,jpj,jpk) , htau  (jpi,jpj)    , dissl(jpi,jpj,jpk) ,     &  
    118          &      avt_k (jpi,jpj,jpk) , avm_k (jpi,jpj,jpk),                          & 
    119          &      avmu_k(jpi,jpj,jpk) , avmv_k(jpi,jpj,jpk), STAT= zdf_tke_alloc      ) 
     118         &      htau  (jpi,jpj)    , dissl(jpi,jpj,jpk) , STAT= zdf_tke_alloc      ) 
    120119         ! 
    121120      IF( lk_mpp             )   CALL mpp_sum ( zdf_tke_alloc ) 
     
    178177         avmu(:,:,:) = avmu_k(:,:,:)  
    179178         avmv(:,:,:) = avmv_k(:,:,:)  
    180       ENDIF  
     179      ENDIF 
     180      ! 
     181      IF( nn_spp_tkelc > 0 ) THEN 
     182          rn_lc0 = rn_lc 
     183          CALL spp_gen(kt,rn_lc0,nn_spp_tkelc,rn_tkelc_sd,    jk_spp_tkelc ) 
     184      ENDIF 
     185      IF( nn_spp_tkedf > 0 ) THEN 
     186          rn_ediff0 = rn_ediff 
     187          CALL spp_gen(kt,rn_ediff0,nn_spp_tkedf,rn_tkedf_sd, jk_spp_tkedf ) 
     188      ENDIF 
     189      IF( nn_spp_tkeds > 0 ) THEN 
     190          rn_ediss0 = rn_ediss 
     191          CALL spp_gen(kt,rn_ediss0,nn_spp_tkeds,rn_tkeds_sd, jk_spp_tkeds ) 
     192      ENDIF 
     193      IF( nn_spp_tkebb > 0 ) THEN 
     194          rn_ebb0 = rn_ebb 
     195          CALL spp_gen(kt,rn_ebb0,nn_spp_tkebb,rn_tkebb_sd,   jk_spp_tkebb ) 
     196      ENDIF 
     197      IF( nn_spp_tkefr > 0 ) THEN 
     198          rn_efr0 = rn_efr 
     199          CALL spp_gen(kt,rn_efr0,nn_spp_tkefr,rn_tkefr_sd,   jk_spp_tkefr ) 
     200      ENDIF 
    181201      ! 
    182202      CALL tke_tke      ! now tke (en) 
     
    188208      avmu_k(:,:,:) = avmu(:,:,:)  
    189209      avmv_k(:,:,:) = avmv(:,:,:)  
     210      ! 
     211      IF (  kt .eq. nitend ) THEN 
     212        DEALLOCATE ( rn_lc0, rn_ediff0, rn_ediss0, rn_ebb0, rn_efr0 ) 
     213      ENDIF 
    190214      ! 
    191215   END SUBROUTINE zdf_tke 
     
    214238      REAL(wp) ::   zrhoa  = 1.22                   ! Air density kg/m3 
    215239      REAL(wp) ::   zcdrag = 1.5e-3                 ! drag coefficient 
    216       REAL(wp) ::   zbbrau, zesh2                   ! temporary scalars 
    217       REAL(wp) ::   zfact1, zfact2, zfact3          !    -         - 
     240      REAL(wp) ::   zesh2                           ! temporary scalars 
     241      REAL(wp) ::   zfact1                          !    -         - 
    218242      REAL(wp) ::   ztx2  , zty2  , zcof            !    -         - 
    219243      REAL(wp) ::   ztau  , zdif                    !    -         - 
     
    222246!!bfr      REAL(wp) ::   zebot                           !    -         - 
    223247      INTEGER , POINTER, DIMENSION(:,:  ) :: imlc 
    224       REAL(wp), POINTER, DIMENSION(:,:  ) :: zhlc 
     248      REAL(wp), POINTER, DIMENSION(:,:  ) :: zhlc, zbbrau,zfact2,zfact3 
    225249      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw 
    226250      !!-------------------------------------------------------------------- 
     
    229253      ! 
    230254      CALL wrk_alloc( jpi,jpj, imlc )    ! integer 
    231       CALL wrk_alloc( jpi,jpj, zhlc )  
     255      CALL wrk_alloc( jpi,jpj, zhlc ) 
     256      CALL wrk_alloc( jpi,jpj, zbbrau ) 
     257      CALL wrk_alloc( jpi,jpj, zfact2 ) 
     258      CALL wrk_alloc( jpi,jpj, zfact3 ) 
    232259      CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw )  
    233260      ! 
    234       zbbrau = rn_ebb / rau0       ! Local constant initialisation 
     261      zbbrau = rn_ebb0 / rau0       ! Local constant initialisation 
    235262      zfact1 = -.5_wp * rdt  
    236       zfact2 = 1.5_wp * rdt * rn_ediss 
    237       zfact3 = 0.5_wp       * rn_ediss 
     263      zfact2 = 1.5_wp * rdt * rn_ediss0 
     264      zfact3 = 0.5_wp       * rn_ediss0 
    238265      ! 
    239266      ! 
     
    250277      DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0) 
    251278         DO ji = fs_2, fs_jpim1   ! vector opt. 
    252             en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
     279            en(ji,jj,1) = MAX( rn_emin0, zbbrau(ji,jj) * taum(ji,jj) ) * tmask(ji,jj,1) 
    253280         END DO 
    254281      END DO 
     
    315342                  !                                           ! vertical velocity due to LC 
    316343                  zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) ) 
    317                   zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) ) 
     344                  zwlc = zind * rn_lc0(ji,jj) * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) ) 
    318345                  !                                           ! TKE Langmuir circulation source term 
    319                   en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     346                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc ) /   & 
     347                     &   zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     348 
    320349               END DO 
    321350            END DO 
     
    360389               zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw) 
    361390               zd_lw(ji,jj,jk) = zzd_lw 
    362                zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk) 
     391               zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2(ji,jj) * dissl(ji,jj,jk) * tmask(ji,jj,jk) 
    363392               ! 
    364393               !                                   ! right hand side in en 
    365394               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zesh2  -   avt(ji,jj,jk) * rn2(ji,jj,jk)    & 
    366                   &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) & 
     395                  &                                 + zfact3(ji,jj) * dissl(ji,jj,jk) * en (ji,jj,jk)  ) & 
    367396                  &                                 * wmask(ji,jj,jk) 
    368397            END DO 
     
    420449            DO jj = 2, jpjm1 
    421450               DO ji = fs_2, fs_jpim1   ! vector opt. 
    422                   en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    423                      &                                 * ( 1._wp - fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     451                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr0(ji,jj) * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
     452                     &                                 * MAX(0.,1._wp - 2.*fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    424453               END DO 
    425454            END DO 
     
    429458            DO ji = fs_2, fs_jpim1   ! vector opt. 
    430459               jk = nmln(ji,jj) 
    431                en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    432                   &                                 * ( 1._wp - fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     460               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr0(ji,jj) * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
     461                  &                                 * MAX(0.,1._wp - 2.*fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    433462            END DO 
    434463         END DO 
     
    445474                  zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean  
    446475                  zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications... 
    447                   en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    448                      &                        * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     476                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau(ji,jj) * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
     477                     &                        * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    449478               END DO 
    450479            END DO 
     
    455484      CALL wrk_dealloc( jpi,jpj, imlc )    ! integer 
    456485      CALL wrk_dealloc( jpi,jpj, zhlc )  
     486      CALL wrk_dealloc( jpi,jpj, zbbrau )  
     487      CALL wrk_dealloc( jpi,jpj, zfact2 )  
     488      CALL wrk_dealloc( jpi,jpj, zfact3 )  
    457489      CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw )  
    458490      ! 
     
    637669            DO ji = fs_2, fs_jpim1   ! vector opt. 
    638670               zsqen = SQRT( en(ji,jj,jk) ) 
    639                zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen 
     671               zav   = rn_ediff0(ji,jj) * zmxlm(ji,jj,jk) * zsqen 
    640672               avm  (ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk) 
    641673               avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 
     
    643675            END DO 
    644676         END DO 
     677         IF(nn_spp_avt > 0 ) CALL spp_gen(1 ,avt(:,:,jk),nn_spp_avt,rn_avt_sd, jk_spp_avt, jk) 
     678         IF(nn_spp_avm > 0 ) CALL spp_gen(1 ,avm(:,:,jk),nn_spp_avm,rn_avm_sd, jk_spp_avm, jk) 
    645679      END DO 
    646680      CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged) 
     
    710744      !!---------------------------------------------------------------------- 
    711745      INTEGER ::   ji, jj, jk   ! dummy loop indices 
    712       INTEGER ::   ios 
     746      INTEGER ::   ios, ierr 
    713747      !! 
    714748      NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin  ,   & 
     
    767801         rn_mxl0 = rmxl_min 
    768802      ENDIF 
     803 
     804      ALLOCATE( rn_lc0   (jpi,jpj) ) ; rn_lc0    = rn_lc 
     805      ALLOCATE( rn_ediff0(jpi,jpj) ) ; rn_ediff0 = rn_ediff 
     806      ALLOCATE( rn_ediss0(jpi,jpj) ) ; rn_ediss0 = rn_ediss 
     807      ALLOCATE( rn_ebb0  (jpi,jpj) ) ; rn_ebb0   = rn_ebb 
     808      ALLOCATE( rn_efr0  (jpi,jpj) ) ; rn_efr0   = rn_efr 
    769809       
    770       IF( nn_etau == 2  )   CALL zdf_mxl( nit000 )      ! Initialization of nmln  
     810      IF( nn_etau == 2 ) THEN 
     811          ierr = zdf_mxl_alloc() 
     812          nmln(:,:) = nlb10           ! Initialization of nmln 
     813      ENDIF 
    771814 
    772815      !                               !* depth of penetration of surface tke 
     
    836879              CALL tke_avn                               ! recompute avt, avm, avmu, avmv and dissl (approximation) 
    837880              ! 
    838               avt_k (:,:,:) = avt (:,:,:) 
    839               avm_k (:,:,:) = avm (:,:,:) 
    840               avmu_k(:,:,:) = avmu(:,:,:) 
    841               avmv_k(:,:,:) = avmv(:,:,:) 
    842               ! 
    843881              DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke( jit )   ;   END DO 
    844882           ENDIF 
    845883        ELSE                                   !* Start from rest 
    846884           en(:,:,:) = rn_emin * tmask(:,:,:) 
    847            DO jk = 1, jpk                           ! set the Kz to the background value 
    848               avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 
    849               avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 
    850               avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk) 
    851               avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk) 
    852            END DO 
    853885        ENDIF 
    854886        ! 
     887        avt_k (:,:,:) = avt (:,:,:) 
     888        avm_k (:,:,:) = avm (:,:,:) 
     889        avmu_k(:,:,:) = avmu(:,:,:) 
     890        avmv_k(:,:,:) = avmv(:,:,:) 
     891         
    855892     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file 
    856893        !                                   ! ------------------- 
Note: See TracChangeset for help on using the changeset viewer.