Changeset 11804 for NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_coare.F90

Timestamp:

2019-10-25T17:25:19+02:00 (5 years ago)

Author:

laurent

Message:

Use of "TURB_FLUXES@…" inside sbcblk.F90, positive latent HF & positive cool-skin temperature now allowed!

File:

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

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

@@ -63 +63 @@
       !!---------------------------------------------------------------------
       !!
-      !!  Cool-Skin scheme according to Fairall et al. 1996, revisited for COARE 3.6 (Fairall et al., 2019)
+      !! Cool-skin parameterization, based on Fairall et al., 1996, revisited for COARE 3.6 (Fairall et al., 2019)
       !!
       !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A.,
@@ -78 +78 @@
       !!     *pSST*       bulk SST (taken at depth gdept_1d(1))          [K]
       !!     *pQlat*      surface latent heat flux                       [K]
-      !!
       !!------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQsw   ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQnsol ! non-solar heat flux to the ocean [W/m^2]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pustar  ! friction velocity, temperature and humidity (u*,t*,q*)
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pSST ! bulk SST [K]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQlat  ! latent heat flux [W/m^2]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQsw   ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQnsol ! non-solar heat flux to the ocean [W/m^2]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pustar  ! friction velocity, temperature and humidity (u*,t*,q*)
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pSST ! bulk SST [K]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQlat  ! latent heat flux [W/m^2]
       !!---------------------------------------------------------------------
       INTEGER  :: ji, jj, jc
@@ -92 +91 @@
          DO ji = 1, jpi
 
-            zQabs  = MIN( -0.1_wp , pQnsol(ji,jj) ) ! first guess, we do not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes$
-            !                                       ! also, we ONLY consider when the viscous layer is loosing heat to the atmosphere, we only deal with cool-skin! => hence the "MIN( -0$
+            zQabs = pQnsol(ji,jj) ! first guess of heat flux absorbed within the viscous sublayer of thicknes delta,
+            !                     !   => we DO not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes zdelta...
 
             zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) )
 
             DO jc = 1, 4 ! because implicit in terms of zdelta...
-               zfr    = MAX( 0.137_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption, Eq.16 (Fairall al. 1996b) /  !LB: why 0.065 and not 0.137 like in the paper??? Beljaars & Zeng use 0.065, not 0.137 !
-               zQabs  = MIN( -0.1_wp , pQnsol(ji,jj) + zfr*pQsw(ji,jj) ) ! Total cooling at the interface
+               zfr = MAX( 0.137_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption, Eq.16 (Fairall al. 1996b) /  !LB: why 0.065 and not 0.137 like in the paper??? Beljaars & Zeng use 0.065, not 0.137 !
+               zQabs = pQnsol(ji,jj) + zfr*pQsw(ji,jj)
                zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) )
             END DO
 
-            dT_cs(ji,jj) = MIN( zQabs*zdelta/rk0_w , 0._wp )   ! temperature increment
+            dT_cs(ji,jj) = zQabs*zdelta/rk0_w   ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!)
 
          END DO
@@ -189 +188 @@
                   &        + 0.45*12.82*(1-EXP(-zHwl/12.82)) ) / zHwl
                zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj) ! first guess of tot. heat flux absorbed in warm layer !LOLO: depends of zfr, which is wild guess... Wrong!!!
-               !PRINT *, '#LBD:  Initial Qsw & Qnsol:', NINT(pQsw(ji,jj)), NINT(pQnsol(ji,jj))
-               !PRINT *, '#LBD:       =>Qabs:', zQabs,' zfr=', zfr
 
                IF ( (ABS(zdTwl) < 1.E-6_wp) .AND. (zQabs <= 0._wp) ) THEN
@@ -196 +193 @@
                   ! since zQabs <= 0._wp
                   ! => no need to go further
-                  !PRINT *, '#LBD: we have not started to to build a WL yet (dT==0)'
-                  !PRINT *, '#LBD: and theres no way it can occur now since zQabs=', zQabs
-                  !PRINT *, '#LBD: => leaving without changing anything...'
                   l_exit = .TRUE.
                END IF
@@ -207 +201 @@
             !LOLO: remove??? has a strong influence !!!
             IF ( (.NOT.(l_exit)) .AND. (Qnt_ac(ji,jj) + zQabs*rdt <= 0._wp) ) THEN
-               !PRINT *, '#LBD: Oh boy! Next Qnt_ac looking weak! =>', Qnt_ac(ji,jj) + zQabs*rdt
-               !PRINT *, '#LBD:  => time to destroy the warm-layer!'
                l_exit       = .TRUE.
                l_destroy_wl = .TRUE.
@@ -219 +211 @@
                ! 1/ A warm layer already exists (dT>0) but it is cooling down because Qabs<0
                ! 2/ Regardless of WL formed (dT==0 or dT>0), we are in the process to initiate one or warm further it !
-
-               !PRINT *, '#LBD:======================================================'
-               !PRINT *, '#LBD: WL action makes sense now! => zQabs,dT_wl=', REAL(zQabs,4), REAL(zdTwl,4)
-               !PRINT *, '#LBD:======================================================'
-               !PRINT *, '#LBD: current values for Qac and Tac=', REAL(Qnt_ac(ji,jj),4), REAL(Tau_ac(ji,jj),4)
 
                ztac = Tau_ac(ji,jj) + MAX(.002_wp , pTau(ji,jj))*rdt      ! updated momentum integral
@@ -238 +225 @@
                   zHwl = MAX( MIN( Hwl_max , zcd1*ztac/SQRT(zqac)) , 0.1_wp ) ! Warm-layer depth
                END DO
-               !PRINT *, '#LBD: updated absorption and WL depth=',  REAL(zfr,4), REAL(zHwl,4)
-               !PRINT *, '#LBD: updated value for Qabs=',  REAL(zQabs,4), 'W/m2'
-               !PRINT *, '#LBD: updated value for Qac =',  REAL(zqac,4), 'J'
 
                IF ( zqac <= 0._wp ) THEN
@@ -263 +247 @@
             END IF
 
-            !PRINT *, '#LBD: exit values for Qac & Tac:', REAL(zqac,4), REAL(ztac,4)
-
             IF ( iwait == 0 ) THEN
                !! Iteration loop within bulk algo is over, time to update what needs to be updated:
                dT_wl(ji,jj)  = zdTwl
                Hz_wl(ji,jj)  = zHwl
-               !PRINT *, '#LBD: FINAL EXIT values for dT_wl & Hz_wl:', REAL(dT_wl(ji,jj),4), REAL(Hz_wl(ji,jj),4)
                Qnt_ac(ji,jj) = zqac ! Updating Qnt_ac, heat integral
                Tau_ac(ji,jj) = ztac
-               !PRINT *, '#LBD: FINAL EXIT values for Qac & Tac:', REAL(Qnt_ac(ji,jj),4), REAL(Tau_ac(ji,jj),4)
-               !PRINT *, '#LBD'
             END IF
 
@@ -284 +263 @@
 
 
-   FUNCTION delta_skin_layer( palpha, pQabs, pQlat, pustar_a )
+   FUNCTION delta_skin_layer( palpha, pQd, pQlat, pustar_a )
       !!---------------------------------------------------------------------
       !! Computes the thickness (m) of the viscous skin layer.
@@ -298 +277 @@
       REAL(wp)                :: delta_skin_layer
       REAL(wp), INTENT(in)    :: palpha   ! thermal expansion coefficient of sea-water (SST accurate enough!)
-      REAL(wp), INTENT(in)    :: pQabs    ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2] => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996
+      REAL(wp), INTENT(in)    :: pQd    ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2] => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996
       REAL(wp), INTENT(in)    :: pQlat    ! latent heat flux [W/m^2]
       REAL(wp), INTENT(in)    :: pustar_a ! friction velocity in the air (u*) [m/s]
       !!---------------------------------------------------------------------
-      REAL(wp) :: zusw, zusw2, zlamb, zQb
-      !!---------------------------------------------------------------------
-
-      zQb = pQabs + 0.026*MIN(pQlat,0._wp)*rCp0_w/rLevap/palpha   ! LOLO: Double check sign + division by palpha !!! units are okay!
-
+      REAL(wp) :: zusw, zusw2, zlamb, zQd, ztf, ztmp
+      !!---------------------------------------------------------------------
+      
+      zQd = pQd + 0.026*MIN(pQlat,0._wp)*rCp0_w/rLevap/palpha   ! LOLO: Double check sign + division by palpha !!! units are okay!
+
+      ztf = 0.5_wp + SIGN(0.5_wp, zQd)  ! Qabs < 0 => cooling of the viscous layer => ztf = 0 (regular case)
+      !                                 ! Qabs > 0 => warming of the viscous layer => ztf = 1 (ex: weak evaporation and strong positive sensible heat flux)
+      !
       zusw  = MAX(pustar_a, 1.E-4_wp) * sq_radrw    ! u* in the water
       zusw2 = zusw*zusw
-
-      zlamb = 6._wp*( 1._wp + (palpha*zcon0/(zusw2*zusw2)*zQb)**0.75 )**(-1./3.) ! see eq.(14) in Fairall et al., 1996
-
-      delta_skin_layer = zlamb*rnu0_w/zusw
-
+      !
+      zlamb = 6._wp*( 1._wp + MAX(palpha*zcon0/(zusw2*zusw2)*zQd, 0._wp)**0.75 )**(-1./3.) ! see Eq.(14) in Fairall et al., 1996
+      !  => zlamb is not used when Qd > 0, and since zcon0 < 0, we just use this "MAX" to prevent FPE errors (something_negative)**0.75
+      !
+      ztmp = rnu0_w/zusw
+      delta_skin_layer = (1._wp-ztf) *     zlamb*ztmp           &  ! regular case, Qd < 0, see Eq.(12) in Fairall et al., 1996
+         &               +   ztf     * MIN(6._wp*ztmp , 0.007_wp)  ! when Qd > 0
    END FUNCTION delta_skin_layer