Skip to content
Snippets Groups Projects
ablmod.F90 63 KiB
Newer Older
Guillaume Samson's avatar
Guillaume Samson committed
            zlup( ji, jpka ) = mxl_min
            zldw( ji, jpka ) = mxl_min
         END DO

         DO jk = 2,jpka-1
            DO ji = 1, jpi
               zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk)
               zldw(ji,jk) = ghw_abl(jk  ) - ghw_abl( 1)
            END DO
         END DO
         !!
         !! BL89 search for lup
         !! ----------------------------------------------------------
         DO jk=2,jpka-1
            !
            DO ji = 1, jpi
               zCF(ji,1:jpka) = 0._wp
               zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a )
               ln_foundl(ji ) = .false.
            END DO
            !
            DO jkup=jk+1,jpka-1
               DO ji = 1, jpi
                  zbuoy1 = MAX( zbn2(ji,jj,jkup  ), rsmall )
                  zbuoy2 = MAX( zbn2(ji,jj,jkup-1), rsmall )
                  zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) * &
                     &               ( zbuoy1*(ghw_abl(jkup  )-ghw_abl(jk)) &
                     &               + zbuoy2*(ghw_abl(jkup-1)-ghw_abl(jk)) )
                  IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN
                     zcff2 = ghw_abl(jkup  ) - ghw_abl(jk)
                     zcff1 = ghw_abl(jkup-1) - ghw_abl(jk)
                     zcff  = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) /   &
                        &    (         zCF(ji,jkup) -         zCF(ji,jkup-1) )
                     zlup(ji,jk) = zcff
                     zlup(ji,jk) = ghw_abl(jkup  ) - ghw_abl(jk)
                     ln_foundl(ji) = .true.
                  END IF
               END DO
            END DO
            !
         END DO
         !!
         !! BL89 search for ldwn
         !! ----------------------------------------------------------
         DO jk=2,jpka-1
            !
            DO ji = 1, jpi
               zCF(ji,1:jpka) = 0._wp
               zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a )
               ln_foundl(ji ) = .false.
            END DO
            !
            DO jkdwn=jk-1,1,-1
               DO ji = 1, jpi
                  zbuoy1 = MAX( zbn2(ji,jj,jkdwn+1), rsmall )
                  zbuoy2 = MAX( zbn2(ji,jj,jkdwn  ), rsmall )
                  zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1)  &
                     &               * ( zbuoy1*(ghw_abl(jk)-ghw_abl(jkdwn+1)) &
                                       + zbuoy2*(ghw_abl(jk)-ghw_abl(jkdwn  )) )
                  IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp  .and. .not. ln_foundl(ji) ) THEN
                     zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1)
                     zcff1 = ghw_abl(jk) - ghw_abl(jkdwn  )
                     zcff  = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) /   &
                        &    (         zCF(ji,jkdwn+1) -         zCF(ji,jkdwn) )
                     zldw(ji,jk) = zcff
                     zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn  )
                     ln_foundl(ji) = .true.
                  END IF
               END DO
            END DO
            !
         END DO

         DO jk = 1, jpka
            DO ji = 1, jpi
               !zcff = 2.*SQRT(2.)*(  zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp)  )**(-3._wp/2._wp)
               zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
               mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min )
               mxld_abl( ji, jj, jk ) = MAX( MIN( zldw( ji, jk ),  zlup( ji, jk ) ), mxl_min )
            END DO
         END DO

      END DO
#   undef zlup
#   undef zldw
         !
     CASE ( 3 )           ! Bougeault & Lacarrere 89 length-scale
         !
#   define zlup zRH
#   define zldw zFC
! zCF is used for matrix inversion
!
       DO jj = 1, jpj      ! outer loop
          !
          DO jk = 2, jpkam1
             DO ji = 1,jpi
                            zcff        = 1.0_wp / e3w_abl( jk )**2
                zdU         = zcff* (u_abl( ji, jj, jk+1, tind)-u_abl( ji, jj, jk  , tind) )**2
                zdV         = zcff* (v_abl( ji, jj, jk+1, tind)-v_abl( ji, jj, jk  , tind) )**2
                zsh2(ji,jk) = SQRT(zdU+zdV)   !<-- zsh2 = SQRT ( ( du/dz )^2 + ( dv/dz )^2 )
                         ENDDO
                  ENDDO
          zsh2(:,      1) = zsh2( :,      2)
          zsh2(:,   jpka) = zsh2( :, jpkam1)

                 DO ji = 1, jpi
            zcff              = zrough(ji,jj) * rn_Lsfc
                        zlup( ji,    1 )  = zcff
            zldw( ji,    1 )  = zcff
            zlup( ji, jpka ) = mxl_min
            zldw( ji, jpka ) = mxl_min
         END DO

         DO jk = 2,jpka-1
            DO ji = 1, jpi
               zlup(ji,jk) = ghw_abl(jpka) - ghw_abl(jk)
               zldw(ji,jk) = ghw_abl(jk  ) - ghw_abl( 1)
            END DO
         END DO
         !!
         !! BL89 search for lup
         !! ----------------------------------------------------------
         DO jk=2,jpka-1
            !
            DO ji = 1, jpi
               zCF(ji,1:jpka) = 0._wp
               zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a )
               ln_foundl(ji ) = .false.
            END DO
            !
            DO jkup=jk+1,jpka-1
               DO ji = 1, jpi
                  zbuoy1             = MAX( zbn2(ji,jj,jkup  ), rsmall )
                  zbuoy2             = MAX( zbn2(ji,jj,jkup-1), rsmall )
                  zCF (ji,jkup) = zCF (ji,jkup-1) + 0.5_wp * e3t_abl(jkup) *                 &
                     &               ( zbuoy1*(ghw_abl(jkup  )-ghw_abl(jk))      &
                     &               + zbuoy2*(ghw_abl(jkup-1)-ghw_abl(jk)) )    &
                                         &                            + 0.5_wp * e3t_abl(jkup) * rn_Rod *        &
                                         &               ( SQRT(tke_abl( ji, jj, jkup  , nt_a ))*zsh2(ji,jkup  ) &
                                         &               + SQRT(tke_abl( ji, jj, jkup-1, nt_a ))*zsh2(ji,jkup-1) )

                  IF( zCF (ji,jkup) * zCF (ji,jkup-1) .le. 0._wp .and. .not. ln_foundl(ji) ) THEN
                     zcff2 = ghw_abl(jkup  ) - ghw_abl(jk)
                     zcff1 = ghw_abl(jkup-1) - ghw_abl(jk)
                     zcff  = ( zcff1 * zCF(ji,jkup) - zcff2 * zCF(ji,jkup-1) ) /   &
                        &    (         zCF(ji,jkup) -         zCF(ji,jkup-1) )
                     zlup(ji,jk) = zcff
                     zlup(ji,jk) = ghw_abl(jkup  ) - ghw_abl(jk)
                     ln_foundl(ji) = .true.
                  END IF
               END DO
            END DO
            !
         END DO
         !!
         !! BL89 search for ldwn
         !! ----------------------------------------------------------
         DO jk=2,jpka-1
            !
            DO ji = 1, jpi
               zCF(ji,1:jpka) = 0._wp
               zCF(ji,  jk  ) = - tke_abl( ji, jj, jk, nt_a )
               ln_foundl(ji ) = .false.
            END DO
            !
            DO jkdwn=jk-1,1,-1
               DO ji = 1, jpi
                  zbuoy1             = MAX( zbn2(ji,jj,jkdwn+1), rsmall )
                  zbuoy2             = MAX( zbn2(ji,jj,jkdwn  ), rsmall )
                  zCF (ji,jkdwn) = zCF (ji,jkdwn+1) + 0.5_wp * e3t_abl(jkdwn+1)  &
                     &               * (zbuoy1*(ghw_abl(jk)-ghw_abl(jkdwn+1))    &
                     &                 +zbuoy2*(ghw_abl(jk)-ghw_abl(jkdwn  )) )  &
                                         &                            + 0.5_wp * e3t_abl(jkup) * rn_Rod *          &
                                         &               ( SQRT(tke_abl( ji, jj, jkdwn+1, nt_a ))*zsh2(ji,jkdwn+1) &
                                         &               + SQRT(tke_abl( ji, jj, jkdwn  , nt_a ))*zsh2(ji,jkdwn  ) )

                  IF(zCF (ji,jkdwn) * zCF (ji,jkdwn+1) .le. 0._wp  .and. .not. ln_foundl(ji) ) THEN
                     zcff2 = ghw_abl(jk) - ghw_abl(jkdwn+1)
                     zcff1 = ghw_abl(jk) - ghw_abl(jkdwn  )
                     zcff  = ( zcff1 * zCF(ji,jkdwn+1) - zcff2 * zCF(ji,jkdwn) ) /   &
                        &    (         zCF(ji,jkdwn+1) -         zCF(ji,jkdwn) )
                     zldw(ji,jk) = zcff
                     zldw(ji,jk) = ghw_abl(jk) - ghw_abl(jkdwn  )
                     ln_foundl(ji) = .true.
                  END IF
               END DO
            END DO
            !
         END DO

         DO jk = 1, jpka
            DO ji = 1, jpi
               !zcff = 2.*SQRT(2.)*(  zldw( ji, jk )**(-2._wp/3._wp) + zlup( ji, jk )**(-2._wp/3._wp)  )**(-3._wp/2._wp)
               zcff = SQRT( zldw( ji, jk ) * zlup( ji, jk ) )
               mxlm_abl( ji, jj, jk ) = MAX( zcff, mxl_min )
               mxld_abl( ji, jj, jk ) = MAX(  MIN( zldw( ji, jk ),  zlup( ji, jk ) ), mxl_min )
            END DO
         END DO

      END DO
#   undef zlup
#   undef zldw
         !
         !
      END SELECT
      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      !                            !  Finalize the computation of turbulent visc./diff.
      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

      !-------------
      DO jj = 1, jpj     ! outer loop
      !-------------
         DO jk = 1, jpka
            DO ji = 1, jpi  ! vector opt.
               zcff  = MAX( rn_phimax, rn_Ric * mxlm_abl( ji, jj, jk ) * mxld_abl( ji, jj, jk )  &
               &     * MAX( zbn2(ji, jj, jk), rsmall ) / tke_abl( ji, jj, jk, nt_a ) )
               zcff2 =  1. / ( 1. + zcff )   !<-- phi_z(z)
               zcff  = mxlm_abl( ji, jj, jk ) * SQRT( tke_abl( ji, jj, jk, nt_a ) )
               !!FL: MAX function probably useless because of the definition of mxl_min
               Avm_abl( ji, jj, jk ) = MAX( rn_Cm * zcff         , avm_bak   )
               Avt_abl( ji, jj, jk ) = MAX( rn_Ct * zcff * zcff2 , avt_bak   )
            END DO
         END DO
      !-------------
      END DO
      !-------------

!---------------------------------------------------------------------------------------------------
   END SUBROUTINE abl_zdf_tke
!===================================================================================================


!===================================================================================================
   SUBROUTINE smooth_pblh( pvar2d, msk )
!---------------------------------------------------------------------------------------------------

      !!----------------------------------------------------------------------
      !!                   ***  ROUTINE smooth_pblh  ***
      !!
      !! ** Purpose :   2D Hanning filter on atmospheric PBL height
      !!
      !! ---------------------------------------------------------------------
      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: msk
      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pvar2d
      INTEGER                                     :: ji,jj
      REAL(wp)                                    :: smth_a, smth_b
      REAL(wp), DIMENSION(jpi,jpj)                :: zdX,zdY,zFX,zFY
      REAL(wp)                                    :: zumsk,zvmsk
      !!
      !!=========================================================
      !!
      !! Hanning filter
      smth_a = 1._wp / 8._wp
      smth_b = 1._wp / 4._wp
      !
      DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls )
         zumsk = msk(ji,jj) * msk(ji+1,jj)
         zdX ( ji, jj ) = ( pvar2d( ji+1,jj ) - pvar2d( ji  ,jj ) ) * zumsk
      END_2D

      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls-1 )
         zvmsk = msk(ji,jj) * msk(ji,jj+1)
         zdY ( ji, jj ) = ( pvar2d( ji, jj+1 ) - pvar2d( ji  ,jj ) ) * zvmsk
      END_2D

      DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls-1 )
         zFY ( ji, jj  ) =   zdY ( ji, jj   )                        &
            & +  smth_a*  ( (zdX ( ji, jj+1 ) - zdX( ji-1, jj+1 ))   &
            &            -  (zdX ( ji, jj   ) - zdX( ji-1, jj   ))  )
      END_2D

      DO_2D( nn_hls, nn_hls-1, nn_hls-1, nn_hls-1 )
         zFX( ji, jj  ) =    zdX( ji, jj   )                         &
           &    + smth_a*(  (zdY( ji+1, jj ) - zdY( ji+1, jj-1))     &
           &             -  (zdY( ji  , jj ) - zdY( ji  , jj-1)) )
      END_2D

      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
         pvar2d( ji  ,jj ) = pvar2d( ji  ,jj )              &
  &         + msk(ji,jj) * smth_b * (                       &
  &                  zFX( ji, jj ) - zFX( ji-1, jj )        &
  &                 +zFY( ji, jj ) - zFY( ji, jj-1 )  )
      END_2D

!---------------------------------------------------------------------------------------------------
   END SUBROUTINE smooth_pblh
!===================================================================================================

!!======================================================================
END MODULE ablmod