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.
sbcblk_skin_ecmwf.F90 in NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC – NEMO

source: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_ecmwf.F90 @ 11637

Last change on this file since 11637 was 11637, checked in by laurent, 5 years ago

LB: preliminary inclusion of "STATION_ASF" test-case!
File size: 11.3 KB
Line 
1MODULE sbcblk_skin_ecmwf
2   !!======================================================================
3   !!                   ***  MODULE  sbcblk_skin_ecmwf  ***
4   !! Computes:
5   !!   * the surface skin temperature (aka SSST) based on the cool-skin/warm-layer
6   !!     scheme used at ECMWF
7   !!    Using formulation/param. of IFS of ECMWF (cycle 40r1)
8   !!   based on IFS doc (avaible online on the ECMWF's website)
9   !!
10   !!   From routine turb_ecmwf maintained and developed in AeroBulk
11   !!            (https://github.com/brodeau/aerobulk)
12   !!
13   !! ** Author: L. Brodeau, September 2019 / AeroBulk (https://github.com/brodeau/aerobulk)
14   !!----------------------------------------------------------------------
15   !! History :  4.0  !  2016-02  (L.Brodeau)   Original code
16   !!----------------------------------------------------------------------
17
18   !!----------------------------------------------------------------------
19   !!  cswl_ecmwf : computes the surface skin temperature (aka SSST)
20   !!               based on the cool-skin/warm-layer scheme used at ECMWF
21   !!----------------------------------------------------------------------
22   USE oce             ! ocean dynamics and tracers
23   USE dom_oce         ! ocean space and time domain
24   USE phycst          ! physical constants
25   USE sbc_oce         ! Surface boundary condition: ocean fields
26
27   USE sbcblk_phy      !LOLO: all thermodynamics functions, rho_air, q_sat, etc...
28
29   USE lib_mpp         ! distribued memory computing library
30   USE in_out_manager  ! I/O manager
31   USE lib_fortran     ! to use key_nosignedzero
32
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC :: CS_ECMWF, WL_ECMWF
38
39   !! Cool-skin related parameters:
40   ! ...
41
42   !! Warm-layer related parameters:
43   REAL(wp), PARAMETER :: rd0  = 3.        !: Depth scale [m] of warm layer, "d" in Eq.11 (Zeng & Beljaars 2005)
44   REAL(wp), PARAMETER :: zRhoCp_w = rho0_w*rCp0_w
45   REAL(wp), PARAMETER :: rNu0 = 1.0       !:  be closer to COARE3p6 ???!LOLO
46   !REAL(wp), PARAMETER :: rNu0 = 0.5       !: Nu (exponent of temperature profile) Eq.11
47   !                                       !: (Zeng & Beljaars 2005) !: set to 0.5 instead of
48   !                                       !: 0.3 to respect a warming of +3 K in calm
49   !                                       !: condition for the insolation peak of +1000W/m^2
50   !!----------------------------------------------------------------------
51CONTAINS
52
53
54   SUBROUTINE CS_ECMWF( pQsw, pQnsol, pustar, pSST,  pdT )
55      !!---------------------------------------------------------------------
56      !!
57      !!  Cool-Skin scheme according to Fairall et al. 1996
58      !! Warm-Layer scheme according to Zeng & Beljaars, 2005 (GRL)
59      !!  " A prognostic scheme of sea surface skin temperature for modeling and data assimilation "
60      !!     as parameterized in IFS Cy45r1 / E.C.M.W.F.
61      !!     ------------------------------------------------------------------
62      !!
63      !!  **   INPUT:
64      !!     *pQsw*       surface net solar radiation into the ocean     [W/m^2] => >= 0 !
65      !!     *pQnsol*     surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
66      !!     *pustar*     friction velocity u*                           [m/s]
67      !!     *pSST*       bulk SST (taken at depth gdept_1d(1))          [K]
68      !!
69      !!  **  INPUT/OUTPUT:
70      !!     *pdT*  : as input =>  previous estimate of dT cool-skin
71      !!              as output =>  new estimate of dT cool-skin
72      !!
73      !!------------------------------------------------------------------
74      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pQsw   ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2]
75      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pQnsol ! non-solar heat flux to the ocean [W/m^2]
76      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pustar  ! friction velocity, temperature and humidity (u*,t*,q*)
77      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)    :: pSST ! bulk SST [K]
78      !!
79      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pdT    ! dT due to cool-skin effect
80      !!---------------------------------------------------------------------
81      INTEGER  ::   ji, jj     ! dummy loop indices
82      REAL(wp) :: zQnet, zQnsol, zlamb, zdelta, zalpha_w, zfr, &
83         & zusw, zusw2
84      !!---------------------------------------------------------------------
85
86      DO jj = 1, jpj
87         DO ji = 1, jpi
88
89            zalpha_w = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!)
90
91            zQnsol = MAX( 1._wp , - pQnsol(ji,jj) ) ! Non-solar heat loss to the atmosphere
92
93            zusw  = MAX(pustar(ji,jj), 1.E-4_wp)*SQRT(roadrw)    ! u* in the water
94            zusw2 = zusw*zusw
95
96            zlamb = 6._wp*( 1._wp + (zQnsol*zalpha_w*rcst_cs/(zusw2*zusw2 ))**0.75 )**(-1./3.)   ! w.r.t COARE 3p6 => seems to ommit absorbed zfr*Qsw (Qnet i.o. Qnsol) and effect of evap...
97            !                                                                                  ! so zlamb not implicit in terms of zdelta (zfr(delta)), so no need to have last guess of delta as in COARE 3.6...
98            zdelta = zlamb*rnu0_w/zusw
99
100            zfr   = MAX( 0.065_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption / Eq. 8.131 / IFS cy40r1, doc, Part IV,
101            zQnet = MAX( 1._wp , zQnsol - zfr*pQsw(ji,jj) ) ! Total cooling at the interface
102
103            !! Update!
104            pdT(ji,jj) =  MIN( - zQnet*zdelta/rk0_w , 0._wp )   ! temperature increment
105
106         END DO
107      END DO
108
109   END SUBROUTINE CS_ECMWF
110
111
112
113   SUBROUTINE WL_ECMWF( pQsw, pQnsol, pustar, pSST, pdT )
114      !!---------------------------------------------------------------------
115      !!
116      !!  Warm-Layer scheme according to Zeng & Beljaars, 2005 (GRL)
117      !!  " A prognostic scheme of sea surface skin temperature for modeling and data assimilation "
118      !!
119      !!    As included in IFS Cy45r1   /  E.C.M.W.F.
120      !!     ------------------------------------------------------------------
121      !!
122      !!  **   INPUT:
123      !!     *pQsw*       surface net solar radiation into the ocean     [W/m^2] => >= 0 !
124      !!     *pQnsol*     surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
125      !!     *pustar*     friction velocity u*                           [m/s]
126      !!     *pSST*       bulk SST  (taken at depth gdept_1d(1))         [K]
127      !!
128      !!   **  INPUT/OUTPUT:
129      !!     *pdT*  : as input =>  previous estimate of dT warm-layer
130      !!             as output =>  new estimate of dT warm-layer
131      !!
132      !!------------------------------------------------------------------
133      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQsw     ! surface net solar radiation into the ocean [W/m^2]     => >= 0 !
134      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQnsol   ! surface net non-solar heat flux into the ocean [W/m^2] => normally < 0 !
135      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pustar   ! friction velocity [m/s]
136      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pSST     ! bulk SST at depth gdept_1d(1) [K]
137      !
138      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pdT    ! dT due to warm-layer effect => difference between "almost surface (right below viscous layer) and depth of bulk SST
139      !
140      INTEGER :: ji,jj
141      !
142      REAL(wp) :: &
143         & zdz,    & !: thickness of the warm-layer [m]
144         & zalpha_w, & !: thermal expansion coefficient of sea-water
145         & ZSRD,   &
146         & dT_wl,   & ! temp. diff. between "almost surface (right below viscous layer) and bottom of WL
147         & zfr,zdL,zdL2, ztmp, &
148         & ZPHI, &
149         & zus_a, zusw, zusw2, &
150         & flg, zQabs, ZL1, ZL2
151      !!---------------------------------------------------------------------
152
153      DO jj = 1, jpj
154         DO ji = 1, jpi
155
156            zdz = rd0 ! first guess for warm-layer depth (and unique..., less advanced than COARE3p6 !)
157
158            ! dT_wl is the difference between "almost surface (right below viscous layer) and bottom of WL (here zdz)
159            ! pdT         "                          "                                    and depth of bulk SST (here gdept_1d(1))!
160            !! => but of course in general the bulk SST is taken shallower than zdz !!! So correction less pronounced!
161            !! => so here since pdT is difference between surface and gdept_1d(1), need to increase fof dT_wl !
162            flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zdz )               ! => 1 when gdept_1d(1)>zdz (pdT(ji,jj) = dT_wl) | 0 when z_s$
163            dT_wl = pdT(ji,jj) / ( flg + (1._wp-flg)*gdept_1d(1)/zdz )
164
165            zalpha_w = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!)
166
167
168            ! *** zfr = Fraction of solar radiation absorbed in warm layer (-)
169            zfr = 1._wp - 0.28_wp*EXP(-71.5_wp*zdz) - 0.27_wp*EXP(-2.8_wp*zdz) - 0.45_wp*EXP(-0.07_wp*zdz)  !: Eq. 8.157
170
171            zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj)       ! tot heat absorbed in warm layer
172
173            zusw  = MAX(pustar(ji,jj), 1.E-4_wp)*SQRT(roadrw)    ! u* in the water
174            zusw2 = zusw*zusw
175
176
177            !! *** 1st rhs term in eq. 8.156 (IFS doc Cy45r1):
178            ZL1 = zQabs / ( zdz * zRhoCp_w * rNu0 ) * (rNu0 + 1._wp)
179
180
181            !! Buoyancy flux and stability parameter (zdl = -z/L) in water
182            ZSRD = zQabs/zRhoCp_w
183            !
184            flg = 0.5_wp + SIGN(0.5_wp, ZSRD)  ! ZSRD > 0. => 1.  / ZSRD < 0. => 0.
185            ztmp = MAX(dT_wl,0._wp)
186            zdl = (flg+1._wp) * ( zusw2 * SQRT(ztmp/(5._wp*zdz*grav*zalpha_w/rNu0)) ) & ! (dT_wl > 0.0 .AND. ZSRD < 0.0)
187               &  +    flg    *  ZSRD                                                                  !   otherwize
188            !
189            zus_a = MAX( pustar(ji,jj), 1.E-4_wp )
190            zdL = zdz*vkarmn*grav/(roadrw)**1.5_wp*zalpha_w*zdL/(zus_a*zus_a*zus_a)
191
192            !! Stability function Phi_t(-z/L) (zdL is -z/L) :
193            flg = 0.5_wp + SIGN(0.5_wp, zdL)  ! zdl > 0. => 1.  / zdl < 0. => 0.
194            zdL2 = zdL*zdL
195            ZPHI =    flg      * ( 1._wp + (5._wp*zdL + 4._wp*zdL2)/(1._wp + 3._wp*zdL + 0.25_wp*zdL2) ) &  ! (zdL > 0) Takaya et al.
196               & + (flg+1._wp) * ( 1._wp/SQRT(1._wp - 16._wp*(-ABS(zdL))) )        ! (zdl < 0) Eq. 8.136
197            !! FOR zdL > 0.0, old relations:
198            !         ZPHI = 1.+5._wp*zdL                                ! Eq. 8.136 (Large et al. 1994)
199            !         ZPHI = 1.+5.0*(zdL+zdL**2)/(1.0+3.0*zdL+zdL**2) ! SHEBA, Grachev et al. 2007
200
201            !! *** 2nd rhs term in eq. 8.156 (IFS doc Cy45r1):
202            ZL2 = - (rNu0 + 1._wp) * vkarmn * zusw / ( zdz * ZPHI )
203
204            ! Forward time / explicit solving of eq. 8.156 (IFS doc Cy45r1): (f_n+1 == pdT(ji,jj) ; f_n == dT_wl)
205            dT_wl = MAX ( dT_wl + rdt*ZL1 + rdt*ZL2*dT_wl , 0._wp )
206
207            ! dT_wl is the difference between "almost surface (right below viscous layer) and bottom of WL (here zdz)
208            !! => but of course in general the bulk SST is taken shallower than zdz !!! So correction less pronounced!
209
210            flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zdz )               ! => 1 when gdept_1d(1)>zdz (pdT(ji,jj) = dT_wl) | 0 when z_s$
211            pdT(ji,jj) = dT_wl * ( flg + (1._wp-flg)*gdept_1d(1)/zdz )
212
213         END DO
214      END DO
215
216   END SUBROUTINE WL_ECMWF
217
218   !!======================================================================
219END MODULE sbcblk_skin_ecmwf
Note: See TracBrowser for help on using the repository browser.