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

Timestamp:

2019-10-23T16:04:12+02:00 (4 years ago)

Author:

laurent

Message:

LB: solid updates+improvements of cool-skin/warm-layer capabilty of COARE and ECMWF bulk algorithms!

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

@@ -44 +44 @@
    PRIVATE
 
-   PUBLIC :: TURB_ECMWF   ! called by sbcblk.F90
+   PUBLIC :: ECMWF_INIT, TURB_ECMWF
 
    !                   !! ECMWF own values for given constants, taken form IFS documentation...
@@ -55 +55 @@
    REAL(wp), PARAMETER ::   alpha_Q = 0.62    !
 
-   INTEGER , PARAMETER ::   nb_itt = 5        ! number of itterations
+   INTEGER , PARAMETER ::   nb_itt = 10             ! number of itterations
 
    !!----------------------------------------------------------------------
 CONTAINS
 
-   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, pdT_cs,           & ! optionals for cool-skin (and warm-layer)
-      &                      pdT_wl )                              ! optionals for warm-layer only
+
+   SUBROUTINE ecmwf_init(l_use_cs, l_use_wl)
+      !!---------------------------------------------------------------------
+      !!                  ***  FUNCTION ecmwf_init  ***
+      !!
+      !! INPUT :
+      !! -------
+      !!    * l_use_cs : use the cool-skin parameterization
+      !!    * l_use_wl : use the warm-layer parameterization
+      !!---------------------------------------------------------------------
+      LOGICAL , INTENT(in) ::   l_use_cs ! use the cool-skin parameterization
+      LOGICAL , INTENT(in) ::   l_use_wl ! use the warm-layer parameterization
+      INTEGER :: ierr
+      !!---------------------------------------------------------------------
+      IF ( l_use_wl ) THEN
+         ierr = 0
+         ALLOCATE ( dT_wl(jpi,jpj), Hz_wl(jpi,jpj), STAT=ierr )
+         !IF( ierr > 0 ) STOP ' ECMWF_INIT => allocation of dT_wl & Hz_wl failed!'
+         dT_wl(:,:)  = 0._wp
+         Hz_wl(:,:)  = rd0 ! (rd0, constant, = 3m is default for Zeng & Beljaars)
+      END IF
+      !!
+      IF ( l_use_cs ) THEN
+         ierr = 0
+         ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr )
+         !IF( ierr > 0 ) STOP ' ECMWF_INIT => allocation of dT_cs failed!'
+         dT_cs(:,:) = -0.25_wp  ! First guess of skin correction
+      END IF
+   END SUBROUTINE ecmwf_init
+
+
+
+   SUBROUTINE turb_ecmwf( kt, 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, pdT_cs,                                & ! optionals for cool-skin (and warm-layer)
+      &                      pdT_wl, pHz_wl )                                                 ! optionals for warm-layer only
       !!----------------------------------------------------------------------
       !!                      ***  ROUTINE  turb_ecmwf  ***
@@ -80 +112 @@
       !! INPUT :
       !! -------
+      !!    *  kt   : current time step (starts at 1)
       !!    *  zt   : height for temperature and spec. hum. of air            [m]
       !!    *  zu   : height for wind speed (usually 10m)                     [m]
@@ -95 +128 @@
       !!
       !!    *  q_s  : SSQ aka saturation specific humidity at temp. T_s       [kg/kg]
-      !!              -> doesn't need to be given a value if skin temp computed (in case l_use_skin=True)
-      !!              -> MUST be given the correct value if not computing skint temp. (in case l_use_skin=False)
+      !!              -> doesn't need to be given a value if skin temp computed (in case l_use_cs=True or l_use_wl=True)
+      !!              -> MUST be given the correct value if not computing skint temp. (in case l_use_cs=False or l_use_wl=False)
       !!
       !! OPTIONAL INPUT:
@@ -103 +136 @@
       !!    *  rad_lw : downwelling longwave radiation at the surface  (>0)   [W/m^2]
       !!    *  slp    : sea-level pressure                                    [Pa]
+      !!
+      !! OPTIONAL OUTPUT:
+      !! ----------------
       !!    * pdT_cs  : SST increment "dT" for cool-skin correction           [K]
       !!    * pdT_wl  : SST increment "dT" for warm-layer correction          [K]
+      !!    * pHz_wl  : thickness of warm-layer                               [m]
       !!
       !! OUTPUT :
@@ -118 +155 @@
       !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
       !!----------------------------------------------------------------------------------
+      INTEGER,  INTENT(in   )                     ::   kt       ! current time step
       REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m]
       REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                             [m]
@@ -139 +177 @@
       REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   slp      !             [Pa]
       REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_cs
-      !
       REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_wl   !             [K]
+      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pHz_wl   !             [m]
       !
       INTEGER :: j_itt
@@ -162 +200 @@
       !!----------------------------------------------------------------------------------
 
+      IF ( kt == nit000 ) CALL ECMWF_INIT(l_use_cs, l_use_wl)
+
       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
@@ -168 +208 @@
       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
+            PRINT *, ' * PROBLEM ('//TRIM(crtnm)//'): you need to provide Qsw, rad_lw & slp to use cool-skin param!'; STOP
+         END IF
       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
+         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
       END IF
 
@@ -322 +356 @@
             !! 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, &
+            CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, 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(:,:)
+            CALL CS_ECMWF( Qsw, ztmp1, u_star, zsst )  ! Qnsol -> ztmp1
+
+            T_s(:,:) = zsst(:,:) + dT_cs(:,:)*tmask(:,:,1)
+            IF( l_use_wl ) T_s(:,:) = T_s(:,:) + dT_wl(:,:)*tmask(:,:,1)
             q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
 
@@ -335 +369 @@
          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, &
+            CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, U_blk, slp, rad_lw, &
                &                   ztmp1, ztmp2)  ! Qnsol -> ztmp1 / Tau -> ztmp2
-            CALL WL_ECMWF( Qsw, ztmp1, u_star, zsst, pdTw )
+            CALL WL_ECMWF( Qsw, ztmp1, u_star, zsst )
             !! Updating T_s and q_s !!!
-            T_s(:,:) = zsst(:,:) + pdTw(:,:)
-            IF( l_use_cs ) T_s(:,:) = T_s(:,:) + pdTc(:,:)
+            T_s(:,:) = zsst(:,:) + dT_wl(:,:)*tmask(:,:,1)
+            IF( l_use_cs ) T_s(:,:) = T_s(:,:) + dT_cs(:,:)*tmask(:,:,1)
             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
@@ -361 +394 @@
       Cen = vkarmn*vkarmn / (log(zu/z0q)*log(zu/z0q))
 
-      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 .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs
+      IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl
+      IF ( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl
 
       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