Changeset 11615 for NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90

Timestamp:

2019-09-30T15:18:21+02:00 (5 years ago)

Author:

laurent

Message:

LB: new cool-skin and warm-layer parameterizations for ECMWF and COARE3.6, COARE3.0 uses same CSWL param as for COARE3.6.

File:

• NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90 (modified) (12 diffs)

NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90

@@ -39 +39 @@
    USE sbc_oce         ! Surface boundary condition: ocean fields
    USE sbcblk_phy      ! all thermodynamics functions, rho_air, q_sat, etc... !LB
-   USE sbcblk_skin     ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB
+   USE sbcblk_skin_ecmwf ! cool-skin/warm layer scheme !LB
 
    IMPLICIT NONE
@@ -60 +60 @@
 CONTAINS
 
-   SUBROUTINE turb_ecmwf( zt, zu, T_s, t_zt, q_s, q_zt, U_zu, &
+   SUBROUTINE turb_ecmwf( zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl, &
       &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,      &
       &                      Cdn, Chn, Cen,                      &
-      &                      Qsw, rad_lw, slp,                   &
-      &                      Tsk_b )
+      &                      Qsw, rad_lw, slp, pdT_cs,           & ! optionals for cool-skin (and warm-layer)
+      &                      pdT_wl )                              ! optionals for warm-layer only
       !!----------------------------------------------------------------------
-       !!                      ***  ROUTINE  turb_ecmwf  ***
+      !!                      ***  ROUTINE  turb_ecmwf  ***
       !!
       !! ** Purpose :   Computes turbulent transfert coefficients of surface
-      !!                fluxes according to IFS doc. (cycle 40)
+      !!                fluxes according to IFS doc. (cycle 45r1)
       !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
       !!                Returns the effective bulk wind speed at zu to be used in the bulk formulas
@@ -82 +82 @@
       !!    *  zt   : height for temperature and spec. hum. of air            [m]
       !!    *  zu   : height for wind speed (usually 10m)                     [m]
-      !!    *  U_zu : scalar wind speed at zu                                 [m/s]
       !!    *  t_zt : potential air temperature at zt                         [K]
       !!    *  q_zt : specific humidity of air at zt                          [kg/kg]
+      !!    *  U_zu : scalar wind speed at zu                                 [m/s]
+      !!    * l_use_cs : use the cool-skin parameterization
+      !!    * l_use_wl : use the warm-layer parameterization
       !!
       !! INPUT/OUTPUT:
@@ -101 +103 @@
       !!    *  rad_lw : downwelling longwave radiation at the surface  (>0)   [W/m^2]
       !!    *  slp    : sea-level pressure                                    [Pa]
-      !!    *  Tsk_b  : estimate of skin temperature at previous time-step    [K]
+      !!    * pdT_cs  : SST increment "dT" for cool-skin correction           [K]
+      !!    * pdT_wl  : SST increment "dT" for warm-layer correction          [K]
       !!
       !! OUTPUT :
@@ -122 +125 @@
       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity at zt             [kg/kg]
       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu                [m/s]
+      LOGICAL , INTENT(in   )                     ::   l_use_cs ! use the cool-skin parameterization
+      LOGICAL , INTENT(in   )                     ::   l_use_wl ! use the warm-layer parameterization
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau)
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
@@ -133 +138 @@
       REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   rad_lw   !             [W/m^2]
       REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   slp      !             [Pa]
-      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   Tsk_b    !             [Pa]
+      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_cs
+      !
+      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_wl   !             [K]
       !
       INTEGER :: j_itt
@@ -145 +152 @@
          &  z0, z0t, z0q
       !
-      ! Cool skin:
-      LOGICAL :: l_use_skin = .FALSE.
-      REAL(wp), DIMENSION(jpi,jpj) :: zsst    ! to back up the initial bulk SST
+      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
+         &                zsst,   &  ! to back up the initial bulk SST
+         &                pdTc,   &  ! SST increment "dT" for cool-skin correction           [K]
+         &                pdTw       ! SST increment "dT" for warm layer correction          [K]
       !
       REAL(wp), DIMENSION(jpi,jpj) ::   func_m, func_h
       REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2
-      !!----------------------------------------------------------------------------------
-
-      ! Cool skin ?
-      IF( PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) ) THEN
-         l_use_skin = .TRUE.
-      END IF
-      IF (lwp) PRINT *, ' *** LOLO(sbcblk_algo_ecmwf.F90) => l_use_skin =', l_use_skin
-
-      ! Identical first gess as in COARE, with IFS parameter values though...
-      !
+      CHARACTER(len=40), PARAMETER :: crtnm = 'turb_ecmwf@sbcblk_algo_ecmwf.F90'
+      !!----------------------------------------------------------------------------------
+
       l_zt_equal_zu = .FALSE.
       IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
 
-      !! Initialization for cool skin:
-      zsst   = T_s    ! save the bulk SST
-      IF( l_use_skin ) THEN
-         ! First guess for skin temperature:
-         IF( PRESENT(Tsk_b) ) THEN
-            T_s = Tsk_b
-         ELSE
-            T_s = T_s - 0.25     ! sst - 0.25
-         END IF
-         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s
+      !! Initializations for cool skin and warm layer:
+      IF ( l_use_cs ) THEN
+         IF( .NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp)) ) THEN
+            PRINT *, ' * PROBLEM ('//trim(crtnm)//'): you need to provide Qsw, rad_lw & slp to use cool-skin param!'
+            STOP
+         END IF
+         ALLOCATE ( pdTc(jpi,jpj) )
+         pdTc(:,:) = -0.25_wp  ! First guess of skin correction
       END IF
 
+      IF ( l_use_wl ) THEN
+         IF(.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) THEN
+            PRINT *, ' * PROBLEM ('//trim(crtnm)//'): you need to provide Qsw, rad_lw & slp to use warm-layer param!'
+            STOP
+         END IF
+         ALLOCATE ( pdTw(jpi,jpj) )
+         pdTw(:,:) = 0._wp
+      END IF
+
+      IF ( l_use_cs .OR. l_use_wl ) THEN
+         ALLOCATE ( zsst(jpi,jpj) )
+         zsst = T_s ! backing up the bulk SST
+         IF( l_use_cs ) T_s = T_s - 0.25_wp   ! First guess of correction
+         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s !LOLO WL too!!!
+      END IF
+
+
+      ! Identical first gess as in COARE, with IFS parameter values though...
+      !
       !! First guess of temperature and humidity at height zu:
       t_zu = MAX( t_zt ,  180._wp )   ! who knows what's given on masked-continental regions...
@@ -197 +215 @@
       z0t    = MIN(ABS(z0t), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO
 
-      ztmp2  = vkarmn/ztmp0
-      Cd     = ztmp2*ztmp2    ! first guess of Cd
+      Cd     = (vkarmn/ztmp0)**2    ! first guess of Cd
 
       ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd
@@ -265 +282 @@
          z0q    = MIN( ABS( alpha_Q*ztmp1                     ) , 0.001_wp)
 
-         !! Update wind at 10m taking into acount convection-related wind gustiness:
-         !! => Chap. 3.2, IFS doc - Cy40r1, Eq.3.17 and Eq.3.18 + Eq.3.8
-         ! Only true when unstable (L<0) => when ztmp0 < 0 => - !!!
-         ztmp2 = ztmp2 * ( MAX(-zi0*Linv/vkarmn , 0._wp))**(2._wp/3._wp) ! => w*^2  (combining Eq. 3.8 and 3.18, hap.3, IFS doc - Cy31r1)
-         !! => equivalent using Beta=1 (gustiness parameter, 1.25 for COARE, also zi0=600 in COARE..)
-         U_blk = MAX( SQRT(U_zu*U_zu + ztmp2) , 0.2_wp )              ! eq.3.17, Chap.3, p.32, IFS doc - Cy31r1
+         !! Update wind at zu with convection-related wind gustiness in unstable conditions (Chap. 3.2, IFS doc - Cy40r1, Eq.3.17 and Eq.3.18 + Eq.3.8)
+         ztmp2 = Beta0*Beta0*ztmp2*(MAX(-zi0*Linv/vkarmn,0._wp))**(2._wp/3._wp) ! square of wind gustiness contribution  (combining Eq. 3.8 and 3.18, hap.3, IFS doc - Cy31r1)
+         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before zi0
+         U_blk = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed
          ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0.
 
@@ -303 +318 @@
          func_h = log(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv)
 
-         !! SKIN related part
-         !! -----------------
-         IF( l_use_skin ) THEN
-            !! compute transfer coefficients at zu : lolo: verifier...
-            Ch = vkarmn*vkarmn/(func_m*func_h)
-            ztmp1 = LOG(zu) - LOG(z0q) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0q*Linv)   ! func_q
-            Ce = vkarmn*vkarmn/(func_m*ztmp1)
-            ! Non-Solar heat flux to the ocean:
-            ztmp1 = U_blk*MAX(rho_air(t_zu, q_zu, slp), 1._wp)     ! rho*U10
-            ztmp2 = T_s*T_s
-            ztmp1 = ztmp1 * ( Ce*L_vap(T_s)*(q_zu - q_s) + Ch*cp_air(q_zu)*(t_zu - T_s) ) & ! Total turb. heat flux
-               &       +      emiss_w*(rad_lw - stefan*ztmp2*ztmp2)                         ! Net longwave flux
-            !!         => "ztmp1" is the net non-solar surface heat flux !
-            !! Updating the values of the skin temperature T_s and q_s :
-            CALL CSWL_ECMWF( Qsw, ztmp1, u_star, zsst, T_s )
-            q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp)  ! 200 -> just to avoid numerics problem on masked regions if silly values are given
-         END IF
-
-         IF( (l_use_skin).OR.(.NOT. l_zt_equal_zu) ) THEN
+
+         IF( l_use_cs ) THEN
+            !! Cool-skin contribution
+
+            CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, &
+               &                   ztmp1, ztmp0,  Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> ztmp0
+
+            CALL CS_ECMWF( Qsw, ztmp1, u_star, zsst, pdTc )  ! Qnsol -> ztmp1
+
+            T_s(:,:) = zsst(:,:) + pdTc(:,:)
+            IF( l_use_wl ) T_s(:,:) = T_s(:,:) + pdTw(:,:)
+            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
+
+         END IF
+
+         IF( l_use_wl ) THEN
+            !! Warm-layer contribution
+            CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, &
+               &                   ztmp1, ztmp2)  ! Qnsol -> ztmp1 / Tau -> ztmp2
+            CALL WL_ECMWF( Qsw, ztmp1, u_star, zsst, pdTw )
+            !! Updating T_s and q_s !!!
+            T_s(:,:) = zsst(:,:) + pdTw(:,:)
+            IF( l_use_cs ) T_s(:,:) = T_s(:,:) + pdTc(:,:)
+            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
+         END IF
+
+
+         IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN
             dt_zu = t_zu - T_s ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
             dq_zu = q_zu - q_s ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
@@ -337 +361 @@
       Cen = vkarmn*vkarmn / (log(zu/z0q)*log(zu/z0q))
 
-   END SUBROUTINE TURB_ECMWF
+      IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = pdTc
+      IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = pdTw
+
+      IF ( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst )
+      IF (          l_use_cs      ) DEALLOCATE ( pdTc )
+      IF (          l_use_wl      ) DEALLOCATE ( pdTw )
+
+   END SUBROUTINE turb_ecmwf
 
 
@@ -433 +464 @@
 
 
-
-
-
    !!======================================================================
 END MODULE sbcblk_algo_ecmwf