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.F90 in NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/SBC – NEMO

source: NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/SBC/sbcblk.F90 @ 13463

Last change on this file since 13463 was 13463, checked in by andmirek, 4 years ago

Ticket #2195:update to trunk 13461

  • Property svn:keywords set to Id
File size: 91.5 KB
Line 
1MODULE sbcblk
2   !!======================================================================
3   !!                       ***  MODULE  sbcblk  ***
4   !! Ocean forcing:  momentum, heat and freshwater flux formulation
5   !!                         Aerodynamic Bulk Formulas
6   !!                        SUCCESSOR OF "sbcblk_core"
7   !!=====================================================================
8   !! History :  1.0  !  2004-08  (U. Schweckendiek)  Original CORE code
9   !!            2.0  !  2005-04  (L. Brodeau, A.M. Treguier)  improved CORE bulk and its user interface
10   !!            3.0  !  2006-06  (G. Madec)  sbc rewritting
11   !!             -   !  2006-12  (L. Brodeau)  Original code for turb_core
12   !!            3.2  !  2009-04  (B. Lemaire)  Introduce iom_put
13   !!            3.3  !  2010-10  (S. Masson)  add diurnal cycle
14   !!            3.4  !  2011-11  (C. Harris)  Fill arrays required by CICE
15   !!            3.7  !  2014-06  (L. Brodeau)  simplification and optimization of CORE bulk
16   !!            4.0  !  2016-06  (L. Brodeau)  sbcblk_core becomes sbcblk and is not restricted to the CORE algorithm anymore
17   !!                 !                        ==> based on AeroBulk (https://github.com/brodeau/aerobulk/)
18   !!            4.0  !  2016-10  (G. Madec)  introduce a sbc_blk_init routine
19   !!            4.0  !  2016-10  (M. Vancoppenolle)  Introduce conduction flux emulator (M. Vancoppenolle)
20   !!            4.0  !  2019-03  (F. Lemarié & G. Samson)  add ABL compatibility (ln_abl=TRUE)
21   !!----------------------------------------------------------------------
22
23   !!----------------------------------------------------------------------
24   !!   sbc_blk_init  : initialisation of the chosen bulk formulation as ocean surface boundary condition
25   !!   sbc_blk       : bulk formulation as ocean surface boundary condition
26   !!   blk_oce_1     : computes pieces of momentum, heat and freshwater fluxes over ocean for ABL model  (ln_abl=TRUE)
27   !!   blk_oce_2     : finalizes momentum, heat and freshwater fluxes computation over ocean after the ABL step  (ln_abl=TRUE)
28   !!             sea-ice case only :
29   !!   blk_ice_1   : provide the air-ice stress
30   !!   blk_ice_2   : provide the heat and mass fluxes at air-ice interface
31   !!   blk_ice_qcn   : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
32   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
33   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag
34   !!----------------------------------------------------------------------
35   USE oce            ! ocean dynamics and tracers
36   USE dom_oce        ! ocean space and time domain
37   USE phycst         ! physical constants
38   USE fldread        ! read input fields
39   USE sbc_oce        ! Surface boundary condition: ocean fields
40   USE cyclone        ! Cyclone 10m wind form trac of cyclone centres
41   USE sbcdcy         ! surface boundary condition: diurnal cycle
42   USE sbcwave , ONLY :   cdn_wave ! wave module
43   USE sbc_ice        ! Surface boundary condition: ice fields
44   USE lib_fortran    ! to use key_nosignedzero
45#if defined key_si3
46   USE ice     , ONLY :   jpl, a_i_b, at_i_b, rn_cnd_s, hfx_err_dif
47   USE icethd_dh      ! for CALL ice_thd_snwblow
48#endif
49   USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009)
50   USE sbcblk_algo_coare3p0 ! => turb_coare3p0 : COAREv3.0 (Fairall et al. 2003)
51   USE sbcblk_algo_coare3p6 ! => turb_coare3p6 : COAREv3.6 (Fairall et al. 2018 + Edson et al. 2013)
52   USE sbcblk_algo_ecmwf    ! => turb_ecmwf    : ECMWF (IFS cycle 45r1)
53   !
54   USE iom            ! I/O manager library
55   USE in_out_manager ! I/O manager
56   USE lib_mpp        ! distribued memory computing library
57   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
58   USE prtctl         ! Print control
59
60   USE sbcblk_phy     ! a catalog of functions for physical/meteorological parameters in the marine boundary layer, rho_air, q_sat, etc...
61
62
63   IMPLICIT NONE
64   PRIVATE
65
66   PUBLIC   sbc_blk_init  ! called in sbcmod
67   PUBLIC   sbc_blk       ! called in sbcmod
68   PUBLIC   blk_oce_1     ! called in sbcabl
69   PUBLIC   blk_oce_2     ! called in sbcabl
70#if defined key_si3
71   PUBLIC   blk_ice_1     ! routine called in icesbc
72   PUBLIC   blk_ice_2     ! routine called in icesbc
73   PUBLIC   blk_ice_qcn   ! routine called in icesbc
74#endif
75
76   INTEGER , PUBLIC, PARAMETER ::   jp_wndi  =  1   ! index of 10m wind velocity (i-component) (m/s)    at T-point
77   INTEGER , PUBLIC, PARAMETER ::   jp_wndj  =  2   ! index of 10m wind velocity (j-component) (m/s)    at T-point
78   INTEGER , PUBLIC, PARAMETER ::   jp_tair  =  3   ! index of 10m air temperature             (Kelvin)
79   INTEGER , PUBLIC, PARAMETER ::   jp_humi  =  4   ! index of specific humidity               ( % )
80   INTEGER , PUBLIC, PARAMETER ::   jp_qsr   =  5   ! index of solar heat                      (W/m2)
81   INTEGER , PUBLIC, PARAMETER ::   jp_qlw   =  6   ! index of Long wave                       (W/m2)
82   INTEGER , PUBLIC, PARAMETER ::   jp_prec  =  7   ! index of total precipitation (rain+snow) (Kg/m2/s)
83   INTEGER , PUBLIC, PARAMETER ::   jp_snow  =  8   ! index of snow (solid prcipitation)       (kg/m2/s)
84   INTEGER , PUBLIC, PARAMETER ::   jp_slp   =  9   ! index of sea level pressure              (Pa)
85   INTEGER , PUBLIC, PARAMETER ::   jp_uoatm = 10   ! index of surface current (i-component)
86   !                                                !          seen by the atmospheric forcing (m/s) at T-point
87   INTEGER , PUBLIC, PARAMETER ::   jp_voatm = 11   ! index of surface current (j-component)
88   !                                                !          seen by the atmospheric forcing (m/s) at T-point
89   INTEGER , PUBLIC, PARAMETER ::   jp_hpgi  = 12   ! index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point
90   INTEGER , PUBLIC, PARAMETER ::   jp_hpgj  = 13   ! index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point
91   INTEGER , PUBLIC, PARAMETER ::   jpfld    = 13   ! maximum number of files to read
92
93   ! Warning: keep this structure allocatable for Agrif...
94   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input atmospheric fields (file informations, fields read)
95
96   !                           !!* Namelist namsbc_blk : bulk parameters
97   LOGICAL  ::   ln_NCAR        ! "NCAR"      algorithm   (Large and Yeager 2008)
98   LOGICAL  ::   ln_COARE_3p0   ! "COARE 3.0" algorithm   (Fairall et al. 2003)
99   LOGICAL  ::   ln_COARE_3p6   ! "COARE 3.6" algorithm   (Edson et al. 2013)
100   LOGICAL  ::   ln_ECMWF       ! "ECMWF"     algorithm   (IFS cycle 45r1)
101   !
102   LOGICAL  ::   ln_Cd_L12      ! ice-atm drag = F( ice concentration )                        (Lupkes et al. JGR2012)
103   LOGICAL  ::   ln_Cd_L15      ! ice-atm drag = F( ice concentration, atmospheric stability ) (Lupkes et al. JGR2015)
104   !
105   LOGICAL  ::   ln_crt_fbk     ! Add surface current feedback to the wind stress computation  (Renault et al. 2020)
106   REAL(wp) ::   rn_stau_a      ! Alpha and Beta coefficients of Renault et al. 2020, eq. 10: Stau = Alpha * Wnd + Beta
107   REAL(wp) ::   rn_stau_b      !
108   !
109   REAL(wp)         ::   rn_pfac   ! multiplication factor for precipitation
110   REAL(wp), PUBLIC ::   rn_efac   ! multiplication factor for evaporation
111   REAL(wp)         ::   rn_zqt    ! z(q,t) : height of humidity and temperature measurements
112   REAL(wp)         ::   rn_zu     ! z(u)   : height of wind measurements
113   !
114   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   Cdn_oce, Chn_oce, Cen_oce  ! neutral coeffs over ocean (L15 bulk scheme and ABL)
115   REAL(wp),         ALLOCATABLE, DIMENSION(:,:) ::   Cd_ice , Ch_ice , Ce_ice   ! transfert coefficients over ice
116   REAL(wp),         ALLOCATABLE, DIMENSION(:,:) ::   t_zu, q_zu                 ! air temp. and spec. hum. at wind speed height (L15 bulk scheme)
117
118   LOGICAL  ::   ln_skin_cs     ! use the cool-skin (only available in ECMWF and COARE algorithms) !LB
119   LOGICAL  ::   ln_skin_wl     ! use the warm-layer parameterization (only available in ECMWF and COARE algorithms) !LB
120   LOGICAL  ::   ln_humi_sph    ! humidity read in files ("sn_humi") is specific humidity [kg/kg] if .true. !LB
121   LOGICAL  ::   ln_humi_dpt    ! humidity read in files ("sn_humi") is dew-point temperature [K] if .true. !LB
122   LOGICAL  ::   ln_humi_rlh    ! humidity read in files ("sn_humi") is relative humidity     [%] if .true. !LB
123   LOGICAL  ::   ln_tpot        !!GS: flag to compute or not potential temperature
124   !
125   INTEGER  ::   nhumi          ! choice of the bulk algorithm
126   !                            ! associated indices:
127   INTEGER, PARAMETER :: np_humi_sph = 1
128   INTEGER, PARAMETER :: np_humi_dpt = 2
129   INTEGER, PARAMETER :: np_humi_rlh = 3
130
131   INTEGER  ::   nblk           ! choice of the bulk algorithm
132   !                            ! associated indices:
133   INTEGER, PARAMETER ::   np_NCAR      = 1   ! "NCAR" algorithm        (Large and Yeager 2008)
134   INTEGER, PARAMETER ::   np_COARE_3p0 = 2   ! "COARE 3.0" algorithm   (Fairall et al. 2003)
135   INTEGER, PARAMETER ::   np_COARE_3p6 = 3   ! "COARE 3.6" algorithm   (Edson et al. 2013)
136   INTEGER, PARAMETER ::   np_ECMWF     = 4   ! "ECMWF" algorithm       (IFS cycle 45r1)
137
138   !! * Substitutions
139#  include "do_loop_substitute.h90"
140   !!----------------------------------------------------------------------
141   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
142   !! $Id$
143   !! Software governed by the CeCILL license (see ./LICENSE)
144   !!----------------------------------------------------------------------
145CONTAINS
146
147   INTEGER FUNCTION sbc_blk_alloc()
148      !!-------------------------------------------------------------------
149      !!             ***  ROUTINE sbc_blk_alloc ***
150      !!-------------------------------------------------------------------
151      ALLOCATE( t_zu(jpi,jpj)   , q_zu(jpi,jpj)   ,                                      &
152         &      Cdn_oce(jpi,jpj), Chn_oce(jpi,jpj), Cen_oce(jpi,jpj),                    &
153         &      Cd_ice (jpi,jpj), Ch_ice (jpi,jpj), Ce_ice (jpi,jpj), STAT=sbc_blk_alloc )
154      !
155      CALL mpp_sum ( 'sbcblk', sbc_blk_alloc )
156      IF( sbc_blk_alloc /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_alloc: failed to allocate arrays' )
157   END FUNCTION sbc_blk_alloc
158
159
160   SUBROUTINE sbc_blk_init
161      !!---------------------------------------------------------------------
162      !!                    ***  ROUTINE sbc_blk_init  ***
163      !!
164      !! ** Purpose :   choose and initialize a bulk formulae formulation
165      !!
166      !! ** Method  :
167      !!
168      !!----------------------------------------------------------------------
169      INTEGER  ::   jfpr                  ! dummy loop indice and argument
170      INTEGER  ::   ios, ierror, ioptio   ! Local integer
171      !!
172      CHARACTER(len=100)            ::   cn_dir                ! Root directory for location of atmospheric forcing files
173      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                 ! array of namelist informations on the fields to read
174      TYPE(FLD_N) ::   sn_wndi, sn_wndj , sn_humi, sn_qsr      ! informations about the fields to be read
175      TYPE(FLD_N) ::   sn_qlw , sn_tair , sn_prec, sn_snow     !       "                        "
176      TYPE(FLD_N) ::   sn_slp , sn_uoatm, sn_voatm             !       "                        "
177      TYPE(FLD_N) ::   sn_hpgi, sn_hpgj                        !       "                        "
178      INTEGER     ::   ipka                                    ! number of levels in the atmospheric variable
179      NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw ,                &   ! input fields
180         &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_uoatm, sn_voatm,     &
181         &                 sn_hpgi, sn_hpgj,                                          &
182         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF,             &   ! bulk algorithm
183         &                 cn_dir , rn_zqt, rn_zu,                                    &
184         &                 rn_pfac, rn_efac, ln_Cd_L12, ln_Cd_L15, ln_tpot,           &
185         &                 ln_crt_fbk, rn_stau_a, rn_stau_b,                          &   ! current feedback
186         &                 ln_skin_cs, ln_skin_wl, ln_humi_sph, ln_humi_dpt, ln_humi_rlh  ! cool-skin / warm-layer !LB
187      !!---------------------------------------------------------------------
188      !
189      !                                      ! allocate sbc_blk_core array
190      IF( sbc_blk_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' )
191      !
192      !                             !** read bulk namelist
193      READ  ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901)
194901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_blk in reference namelist' )
195      !
196      READ  ( numnam_cfg, namsbc_blk, IOSTAT = ios, ERR = 902 )
197902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_blk in configuration namelist' )
198      !
199      IF(lwm) WRITE( numond, namsbc_blk )
200      !
201      !                             !** initialization of the chosen bulk formulae (+ check)
202      !                                   !* select the bulk chosen in the namelist and check the choice
203      ioptio = 0
204      IF( ln_NCAR      ) THEN
205         nblk =  np_NCAR        ;   ioptio = ioptio + 1
206      ENDIF
207      IF( ln_COARE_3p0 ) THEN
208         nblk =  np_COARE_3p0   ;   ioptio = ioptio + 1
209      ENDIF
210      IF( ln_COARE_3p6 ) THEN
211         nblk =  np_COARE_3p6   ;   ioptio = ioptio + 1
212      ENDIF
213      IF( ln_ECMWF     ) THEN
214         nblk =  np_ECMWF       ;   ioptio = ioptio + 1
215      ENDIF
216      IF( ioptio /= 1 )   CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' )
217
218      !                             !** initialization of the cool-skin / warm-layer parametrization
219      IF( ln_skin_cs .OR. ln_skin_wl ) THEN
220         !! Some namelist sanity tests:
221         IF( ln_NCAR )      &
222            & CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm' )
223         IF( nn_fsbc /= 1 ) &
224            & CALL ctl_stop( 'sbc_blk_init: Please set "nn_fsbc" to 1 when using cool-skin/warm-layer param.')
225      END IF
226
227      IF( ln_skin_wl ) THEN
228         !! Check if the frequency of downwelling solar flux input makes sense and if ln_dm2dc=T if it is daily!
229         IF( (sn_qsr%freqh  < 0.).OR.(sn_qsr%freqh  > 24.) ) &
230            & CALL ctl_stop( 'sbc_blk_init: Warm-layer param. (ln_skin_wl) not compatible with freq. of solar flux > daily' )
231         IF( (sn_qsr%freqh == 24.).AND.(.NOT. ln_dm2dc) ) &
232            & CALL ctl_stop( 'sbc_blk_init: Please set ln_dm2dc=T for warm-layer param. (ln_skin_wl) to work properly' )
233      END IF
234
235      ioptio = 0
236      IF( ln_humi_sph ) THEN
237         nhumi =  np_humi_sph    ;   ioptio = ioptio + 1
238      ENDIF
239      IF( ln_humi_dpt ) THEN
240         nhumi =  np_humi_dpt    ;   ioptio = ioptio + 1
241      ENDIF
242      IF( ln_humi_rlh ) THEN
243         nhumi =  np_humi_rlh    ;   ioptio = ioptio + 1
244      ENDIF
245      IF( ioptio /= 1 )   CALL ctl_stop( 'sbc_blk_init: Choose one and only one type of air humidity' )
246      !
247      IF( ln_dm2dc ) THEN                 !* check: diurnal cycle on Qsr
248         IF( sn_qsr%freqh /= 24. )   CALL ctl_stop( 'sbc_blk_init: ln_dm2dc=T only with daily short-wave input' )
249         IF( sn_qsr%ln_tint ) THEN
250            CALL ctl_warn( 'sbc_blk_init: ln_dm2dc=T daily qsr time interpolation done by sbcdcy module',   &
251               &           '              ==> We force time interpolation = .false. for qsr' )
252            sn_qsr%ln_tint = .false.
253         ENDIF
254      ENDIF
255      !                                   !* set the bulk structure
256      !                                      !- store namelist information in an array
257      !
258      slf_i(jp_wndi ) = sn_wndi    ;   slf_i(jp_wndj ) = sn_wndj
259      slf_i(jp_qsr  ) = sn_qsr     ;   slf_i(jp_qlw  ) = sn_qlw
260      slf_i(jp_tair ) = sn_tair    ;   slf_i(jp_humi ) = sn_humi
261      slf_i(jp_prec ) = sn_prec    ;   slf_i(jp_snow ) = sn_snow
262      slf_i(jp_slp  ) = sn_slp
263      slf_i(jp_uoatm) = sn_uoatm   ;   slf_i(jp_voatm) = sn_voatm
264      slf_i(jp_hpgi ) = sn_hpgi    ;   slf_i(jp_hpgj ) = sn_hpgj
265      !
266      IF( .NOT. ln_abl ) THEN   ! force to not use jp_hpgi and jp_hpgj, should already be done in namelist_* but we never know...
267         slf_i(jp_hpgi)%clname = 'NOT USED'
268         slf_i(jp_hpgj)%clname = 'NOT USED'
269      ENDIF
270      !
271      !                                      !- allocate the bulk structure
272      ALLOCATE( sf(jpfld), STAT=ierror )
273      IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_init: unable to allocate sf structure' )
274      !
275      !                                      !- fill the bulk structure with namelist informations
276      CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' )
277      !
278      DO jfpr= 1, jpfld
279         !
280         IF(   ln_abl    .AND.                                                      &
281            &    ( jfpr == jp_wndi .OR. jfpr == jp_wndj .OR. jfpr == jp_humi .OR.   &
282            &      jfpr == jp_hpgi .OR. jfpr == jp_hpgj .OR. jfpr == jp_tair     )  ) THEN
283            ipka = jpka   ! ABL: some fields are 3D input
284         ELSE
285            ipka = 1
286         ENDIF
287         !
288         ALLOCATE( sf(jfpr)%fnow(jpi,jpj,ipka) )
289         !
290         IF( TRIM(sf(jfpr)%clrootname) == 'NOT USED' ) THEN    !--  not used field  --!   (only now allocated and set to default)
291            IF(     jfpr == jp_slp  ) THEN
292               sf(jfpr)%fnow(:,:,1:ipka) = 101325._wp   ! use standard pressure in Pa
293            ELSEIF( jfpr == jp_prec .OR. jfpr == jp_snow .OR. jfpr == jp_uoatm .OR. jfpr == jp_voatm ) THEN
294               sf(jfpr)%fnow(:,:,1:ipka) = 0._wp        ! no precip or no snow or no surface currents
295            ELSEIF( ( jfpr == jp_hpgi .OR. jfpr == jp_hpgj ) .AND. .NOT. ln_abl ) THEN
296               DEALLOCATE( sf(jfpr)%fnow )              ! deallocate as not used in this case
297            ELSE
298               WRITE(ctmp1,*) 'sbc_blk_init: no default value defined for field number', jfpr
299               CALL ctl_stop( ctmp1 )
300            ENDIF
301         ELSE                                                  !-- used field  --!
302            IF( sf(jfpr)%ln_tint )   ALLOCATE( sf(jfpr)%fdta(jpi,jpj,ipka,2) )   ! allocate array for temporal interpolation
303            !
304            IF( sf(jfpr)%freqh > 0. .AND. MOD( NINT(3600. * sf(jfpr)%freqh), nn_fsbc * NINT(rn_Dt) ) /= 0 )   &
305               &  CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rn_Dt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.',   &
306               &                 '               This is not ideal. You should consider changing either rn_Dt or nn_fsbc value...' )
307         ENDIF
308      END DO
309      !
310      IF( ln_wave ) THEN
311         !Activated wave module but neither drag nor stokes drift activated
312         IF( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) )   THEN
313            CALL ctl_stop( 'STOP',  'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' )
314            !drag coefficient read from wave model definable only with mfs bulk formulae and core
315         ELSEIF(ln_cdgw .AND. .NOT. ln_NCAR )       THEN
316            CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae')
317         ELSEIF(ln_stcor .AND. .NOT. ln_sdw)                             THEN
318            CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')
319         ENDIF
320      ELSE
321         IF( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor )                &
322            &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    &
323            &                  'with drag coefficient (ln_cdgw =T) '  ,                        &
324            &                  'or Stokes Drift (ln_sdw=T) ' ,                                 &
325            &                  'or ocean stress modification due to waves (ln_tauwoc=T) ',      &
326            &                  'or Stokes-Coriolis term (ln_stcori=T)'  )
327      ENDIF
328      !
329      IF( ln_abl ) THEN       ! ABL: read 3D fields for wind, temperature, humidity and pressure gradient
330         rn_zqt = ght_abl(2)          ! set the bulk altitude to ABL first level
331         rn_zu  = ght_abl(2)
332         IF(lwp) WRITE(numout,*)
333         IF(lwp) WRITE(numout,*) '   ABL formulation: overwrite rn_zqt & rn_zu with ABL first level altitude'
334      ENDIF
335      !
336      ! set transfer coefficients to default sea-ice values
337      Cd_ice(:,:) = rCd_ice
338      Ch_ice(:,:) = rCd_ice
339      Ce_ice(:,:) = rCd_ice
340      !
341      IF(lwp) THEN                     !** Control print
342         !
343         WRITE(numout,*)                  !* namelist
344         WRITE(numout,*) '   Namelist namsbc_blk (other than data information):'
345         WRITE(numout,*) '      "NCAR"      algorithm   (Large and Yeager 2008)     ln_NCAR      = ', ln_NCAR
346         WRITE(numout,*) '      "COARE 3.0" algorithm   (Fairall et al. 2003)       ln_COARE_3p0 = ', ln_COARE_3p0
347         WRITE(numout,*) '      "COARE 3.6" algorithm (Fairall 2018 + Edson al 2013)ln_COARE_3p6 = ', ln_COARE_3p6
348         WRITE(numout,*) '      "ECMWF"     algorithm   (IFS cycle 45r1)            ln_ECMWF     = ', ln_ECMWF
349         WRITE(numout,*) '      Air temperature and humidity reference height (m)   rn_zqt       = ', rn_zqt
350         WRITE(numout,*) '      Wind vector reference height (m)                    rn_zu        = ', rn_zu
351         WRITE(numout,*) '      factor applied on precipitation (total & snow)      rn_pfac      = ', rn_pfac
352         WRITE(numout,*) '      factor applied on evaporation                       rn_efac      = ', rn_efac
353         WRITE(numout,*) '         (form absolute (=0) to relative winds(=1))'
354         WRITE(numout,*) '      use ice-atm drag from Lupkes2012                    ln_Cd_L12    = ', ln_Cd_L12
355         WRITE(numout,*) '      use ice-atm drag from Lupkes2015                    ln_Cd_L15    = ', ln_Cd_L15
356         WRITE(numout,*) '      use surface current feedback on wind stress         ln_crt_fbk   = ', ln_crt_fbk
357         IF(ln_crt_fbk) THEN
358         WRITE(numout,*) '         Renault et al. 2020, eq. 10: Stau = Alpha * Wnd + Beta'
359         WRITE(numout,*) '            Alpha                                         rn_stau_a    = ', rn_stau_a
360         WRITE(numout,*) '            Beta                                          rn_stau_b    = ', rn_stau_b
361         ENDIF
362         !
363         WRITE(numout,*)
364         SELECT CASE( nblk )              !* Print the choice of bulk algorithm
365         CASE( np_NCAR      )   ;   WRITE(numout,*) '   ==>>>   "NCAR" algorithm        (Large and Yeager 2008)'
366         CASE( np_COARE_3p0 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.0" algorithm   (Fairall et al. 2003)'
367         CASE( np_COARE_3p6 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.6" algorithm (Fairall 2018+Edson et al. 2013)'
368         CASE( np_ECMWF     )   ;   WRITE(numout,*) '   ==>>>   "ECMWF" algorithm       (IFS cycle 45r1)'
369         END SELECT
370         !
371         WRITE(numout,*)
372         WRITE(numout,*) '      use cool-skin  parameterization (SSST)  ln_skin_cs  = ', ln_skin_cs
373         WRITE(numout,*) '      use warm-layer parameterization (SSST)  ln_skin_wl  = ', ln_skin_wl
374         !
375         WRITE(numout,*)
376         SELECT CASE( nhumi )              !* Print the choice of air humidity
377         CASE( np_humi_sph )   ;   WRITE(numout,*) '   ==>>>   air humidity is SPECIFIC HUMIDITY     [kg/kg]'
378         CASE( np_humi_dpt )   ;   WRITE(numout,*) '   ==>>>   air humidity is DEW-POINT TEMPERATURE [K]'
379         CASE( np_humi_rlh )   ;   WRITE(numout,*) '   ==>>>   air humidity is RELATIVE HUMIDITY     [%]'
380         END SELECT
381         !
382      ENDIF
383
384      IF(lxios_blkw) THEN
385         CALL iom_blk_wgt_init( cbwxios_context )
386         CALL iom_swap( cxios_context )
387      ENDIF
388     
389      !
390   END SUBROUTINE sbc_blk_init
391
392
393   SUBROUTINE sbc_blk( kt )
394      !!---------------------------------------------------------------------
395      !!                    ***  ROUTINE sbc_blk  ***
396      !!
397      !! ** Purpose :   provide at each time step the surface ocean fluxes
398      !!              (momentum, heat, freshwater and runoff)
399      !!
400      !! ** Method  :
401      !!              (1) READ each fluxes in NetCDF files:
402      !!      the wind velocity (i-component) at z=rn_zu  (m/s) at T-point
403      !!      the wind velocity (j-component) at z=rn_zu  (m/s) at T-point
404      !!      the specific humidity           at z=rn_zqt (kg/kg)
405      !!      the air temperature             at z=rn_zqt (Kelvin)
406      !!      the solar heat                              (W/m2)
407      !!      the Long wave                               (W/m2)
408      !!      the total precipitation (rain+snow)         (Kg/m2/s)
409      !!      the snow (solid precipitation)              (kg/m2/s)
410      !!      ABL dynamical forcing (i/j-components of either hpg or geostrophic winds)
411      !!              (2) CALL blk_oce_1 and blk_oce_2
412      !!
413      !!      C A U T I O N : never mask the surface stress fields
414      !!                      the stress is assumed to be in the (i,j) mesh referential
415      !!
416      !! ** Action  :   defined at each time-step at the air-sea interface
417      !!              - utau, vtau  i- and j-component of the wind stress
418      !!              - taum        wind stress module at T-point
419      !!              - wndm        wind speed  module at T-point over free ocean or leads in presence of sea-ice
420      !!              - qns, qsr    non-solar and solar heat fluxes
421      !!              - emp         upward mass flux (evapo. - precip.)
422      !!              - sfx         salt flux due to freezing/melting (non-zero only if ice is present)
423      !!
424      !! ** References :   Large & Yeager, 2004 / Large & Yeager, 2008
425      !!                   Brodeau et al. Ocean Modelling 2010
426      !!----------------------------------------------------------------------
427      INTEGER, INTENT(in) ::   kt   ! ocean time step
428      !!----------------------------------------------------------------------
429      REAL(wp), DIMENSION(jpi,jpj) ::   zssq, zcd_du, zsen, zevp
430      REAL(wp) :: ztmp
431      !!----------------------------------------------------------------------
432      !
433      CALL fld_read( kt, nn_fsbc, sf )             ! input fields provided at the current time-step
434
435      ! Sanity/consistence test on humidity at first time step to detect potential screw-up:
436      IF( kt == nit000 ) THEN
437         IF(lwp) WRITE(numout,*) ''
438#if defined key_agrif
439         IF(lwp) WRITE(numout,*) ' === AGRIF => Sanity/consistence test on air humidity SKIPPED! :( ==='
440#else
441         ztmp = SUM(tmask(:,:,1)) ! number of ocean points on local proc domain
442         IF( ztmp > 8._wp ) THEN ! test only on proc domains with at least 8 ocean points!
443            ztmp = SUM(sf(jp_humi)%fnow(:,:,1)*tmask(:,:,1))/ztmp ! mean humidity over ocean on proc
444            SELECT CASE( nhumi )
445            CASE( np_humi_sph ) ! specific humidity => expect: 0. <= something < 0.065 [kg/kg] (0.061 is saturation at 45degC !!!)
446               IF(  (ztmp < 0._wp) .OR. (ztmp > 0.065)  ) ztmp = -1._wp
447            CASE( np_humi_dpt ) ! dew-point temperature => expect: 110. <= something < 320. [K]
448               IF( (ztmp < 110._wp).OR.(ztmp > 320._wp) ) ztmp = -1._wp
449            CASE( np_humi_rlh ) ! relative humidity => expect: 0. <= something < 100. [%]
450               IF(  (ztmp < 0._wp) .OR.(ztmp > 100._wp) ) ztmp = -1._wp
451            END SELECT
452            IF(ztmp < 0._wp) THEN
453               IF (lwp) WRITE(numout,'("   Mean humidity value found on proc #",i6.6," is: ",f10.5)') narea, ztmp
454               CALL ctl_stop( 'STOP', 'Something is wrong with air humidity!!!', &
455                  &   ' ==> check the unit in your input files'       , &
456                  &   ' ==> check consistence of namelist choice: specific? relative? dew-point?', &
457                  &   ' ==> ln_humi_sph -> [kg/kg] | ln_humi_rlh -> [%] | ln_humi_dpt -> [K] !!!' )
458            END IF
459         END IF
460         IF(lwp) WRITE(numout,*) ' === Sanity/consistence test on air humidity sucessfuly passed! ==='
461#endif
462         IF(lwp) WRITE(numout,*) ''
463      END IF !IF( kt == nit000 )
464      !                                            ! compute the surface ocean fluxes using bulk formulea
465      IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
466         CALL blk_oce_1( kt, sf(jp_wndi )%fnow(:,:,1), sf(jp_wndj )%fnow(:,:,1),   &   !   <<= in
467            &                sf(jp_tair )%fnow(:,:,1), sf(jp_humi )%fnow(:,:,1),   &   !   <<= in
468            &                sf(jp_slp  )%fnow(:,:,1), sst_m, ssu_m, ssv_m,        &   !   <<= in
469            &                sf(jp_uoatm)%fnow(:,:,1), sf(jp_voatm)%fnow(:,:,1),   &   !   <<= in
470            &                sf(jp_qsr  )%fnow(:,:,1), sf(jp_qlw  )%fnow(:,:,1),   &   !   <<= in (wl/cs)
471            &                tsk_m, zssq, zcd_du, zsen, zevp )                         !   =>> out
472
473         CALL blk_oce_2(     sf(jp_tair )%fnow(:,:,1), sf(jp_qsr  )%fnow(:,:,1),   &   !   <<= in
474            &                sf(jp_qlw  )%fnow(:,:,1), sf(jp_prec )%fnow(:,:,1),   &   !   <<= in
475            &                sf(jp_snow )%fnow(:,:,1), tsk_m,                      &   !   <<= in
476            &                zsen, zevp )                                              !   <=> in out
477      ENDIF
478      !
479#if defined key_cice
480      IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN
481         qlw_ice(:,:,1)   = sf(jp_qlw )%fnow(:,:,1)
482         IF( ln_dm2dc ) THEN
483            qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )
484         ELSE
485            qsr_ice(:,:,1) =          sf(jp_qsr)%fnow(:,:,1)
486         ENDIF
487         tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)
488
489         SELECT CASE( nhumi )
490         CASE( np_humi_sph )
491            qatm_ice(:,:) =           sf(jp_humi)%fnow(:,:,1)
492         CASE( np_humi_dpt )
493            qatm_ice(:,:) = q_sat(    sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
494         CASE( np_humi_rlh )
495            qatm_ice(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) !LB: 0.01 => RH is % percent in file
496         END SELECT
497
498         tprecip(:,:)     = sf(jp_prec)%fnow(:,:,1) * rn_pfac
499         sprecip(:,:)     = sf(jp_snow)%fnow(:,:,1) * rn_pfac
500         wndi_ice(:,:)    = sf(jp_wndi)%fnow(:,:,1)
501         wndj_ice(:,:)    = sf(jp_wndj)%fnow(:,:,1)
502      ENDIF
503#endif
504      !
505   END SUBROUTINE sbc_blk
506
507
508   SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, phumi,         &  ! inp
509      &                      pslp , pst  , pu   , pv,            &  ! inp
510      &                      puatm, pvatm, pqsr , pqlw ,         &  ! inp
511      &                      ptsk , pssq , pcd_du, psen, pevp   )   ! out
512      !!---------------------------------------------------------------------
513      !!                     ***  ROUTINE blk_oce_1  ***
514      !!
515      !! ** Purpose :   if ln_blk=T, computes surface momentum, heat and freshwater fluxes
516      !!                if ln_abl=T, computes Cd x |U|, Ch x |U|, Ce x |U| for ABL integration
517      !!
518      !! ** Method  :   bulk formulae using atmospheric fields from :
519      !!                if ln_blk=T, atmospheric fields read in sbc_read
520      !!                if ln_abl=T, the ABL model at previous time-step
521      !!
522      !! ** Outputs : - pssq    : surface humidity used to compute latent heat flux (kg/kg)
523      !!              - pcd_du  : Cd x |dU| at T-points  (m/s)
524      !!              - psen    : Ch x |dU| at T-points  (m/s)
525      !!              - pevp    : Ce x |dU| at T-points  (m/s)
526      !!---------------------------------------------------------------------
527      INTEGER , INTENT(in   )                 ::   kt     ! time step index
528      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pwndi  ! atmospheric wind at U-point              [m/s]
529      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pwndj  ! atmospheric wind at V-point              [m/s]
530      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   phumi  ! specific humidity at T-points            [kg/kg]
531      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   ptair  ! potential temperature at T-points        [Kelvin]
532      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pslp   ! sea-level pressure                       [Pa]
533      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pst    ! surface temperature                      [Celsius]
534      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pu     ! surface current at U-point (i-component) [m/s]
535      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pv     ! surface current at V-point (j-component) [m/s]
536      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   puatm  ! surface current seen by the atm at T-point (i-component) [m/s]
537      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pvatm  ! surface current seen by the atm at T-point (j-component) [m/s]
538      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pqsr   !
539      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pqlw   !
540      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   ptsk   ! skin temp. (or SST if CS & WL not used)  [Celsius]
541      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   pssq   ! specific humidity at pst                 [kg/kg]
542      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   pcd_du ! Cd x |dU| at T-points                    [m/s]
543      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   psen   ! Ch x |dU| at T-points                    [m/s]
544      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   pevp   ! Ce x |dU| at T-points                    [m/s]
545      !
546      INTEGER  ::   ji, jj               ! dummy loop indices
547      REAL(wp) ::   zztmp                ! local variable
548      REAL(wp) ::   zstmax, zstau
549#if defined key_cyclone
550      REAL(wp), DIMENSION(jpi,jpj) ::   zwnd_i, zwnd_j    ! wind speed components at T-point
551#endif
552      REAL(wp), DIMENSION(jpi,jpj) ::   ztau_i, ztau_j    ! wind stress components at T-point
553      REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s]
554      REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K]
555      REAL(wp), DIMENSION(jpi,jpj) ::   zqair             ! specific humidity     of air at z=rn_zqt [kg/kg]
556      REAL(wp), DIMENSION(jpi,jpj) ::   zcd_oce           ! momentum transfert coefficient over ocean
557      REAL(wp), DIMENSION(jpi,jpj) ::   zch_oce           ! sensible heat transfert coefficient over ocean
558      REAL(wp), DIMENSION(jpi,jpj) ::   zce_oce           ! latent   heat transfert coefficient over ocean
559      REAL(wp), DIMENSION(jpi,jpj) ::   zqla              ! latent heat flux
560      REAL(wp), DIMENSION(jpi,jpj) ::   zztmp1, zztmp2
561      !!---------------------------------------------------------------------
562      !
563      ! local scalars ( place there for vector optimisation purposes)
564      !                           ! Temporary conversion from Celcius to Kelvin (and set minimum value far above 0 K)
565      ptsk(:,:) = pst(:,:) + rt0  ! by default: skin temperature = "bulk SST" (will remain this way if NCAR algorithm used!)
566
567      ! ----------------------------------------------------------------------------- !
568      !      0   Wind components and module at T-point relative to the moving ocean   !
569      ! ----------------------------------------------------------------------------- !
570
571      ! ... components ( U10m - U_oce ) at T-point (unmasked)
572#if defined key_cyclone
573      zwnd_i(:,:) = 0._wp
574      zwnd_j(:,:) = 0._wp
575      CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012)
576      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
577         zwnd_i(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj)
578         zwnd_j(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj)
579         ! ... scalar wind at T-point (not masked)
580         wndm(ji,jj) = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) + zwnd_j(ji,jj) * zwnd_j(ji,jj) )
581      END_2D
582#else
583      ! ... scalar wind module at T-point (not masked)
584      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
585         wndm(ji,jj) = SQRT(  pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj)  )
586      END_2D
587#endif
588      ! ----------------------------------------------------------------------------- !
589      !      I   Solar FLUX                                                           !
590      ! ----------------------------------------------------------------------------- !
591
592      ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle                    ! Short Wave
593      zztmp = 1. - albo
594      IF( ln_dm2dc ) THEN
595         qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)
596      ELSE
597         qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1)
598      ENDIF
599
600
601      ! ----------------------------------------------------------------------------- !
602      !     II   Turbulent FLUXES                                                     !
603      ! ----------------------------------------------------------------------------- !
604
605      ! specific humidity at SST
606      pssq(:,:) = rdct_qsat_salt * q_sat( ptsk(:,:), pslp(:,:) )
607
608      IF( ln_skin_cs .OR. ln_skin_wl ) THEN
609         !! Backup "bulk SST" and associated spec. hum.
610         zztmp1(:,:) = ptsk(:,:)
611         zztmp2(:,:) = pssq(:,:)
612      ENDIF
613
614      ! specific humidity of air at "rn_zqt" m above the sea
615      SELECT CASE( nhumi )
616      CASE( np_humi_sph )
617         zqair(:,:) = phumi(:,:)      ! what we read in file is already a spec. humidity!
618      CASE( np_humi_dpt )
619         !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of d_air and slp !' !LBrm
620         zqair(:,:) = q_sat( phumi(:,:), pslp(:,:) )
621      CASE( np_humi_rlh )
622         !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of RH, t_air and slp !' !LBrm
623         zqair(:,:) = q_air_rh( 0.01_wp*phumi(:,:), ptair(:,:), pslp(:,:) ) !LB: 0.01 => RH is % percent in file
624      END SELECT
625      !
626      ! potential temperature of air at "rn_zqt" m above the sea
627      IF( ln_abl ) THEN
628         ztpot = ptair(:,:)
629      ELSE
630         ! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate
631         !    (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2
632         !    (since reanalysis products provide T at z, not theta !)
633         !#LB: because AGRIF hates functions that return something else than a scalar, need to
634         !     use scalar version of gamma_moist() ...
635         IF( ln_tpot ) THEN
636            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
637               ztpot(ji,jj) = ptair(ji,jj) + gamma_moist( ptair(ji,jj), zqair(ji,jj) ) * rn_zqt
638            END_2D
639         ELSE
640            ztpot = ptair(:,:)
641         ENDIF
642      ENDIF
643
644      !! Time to call the user-selected bulk parameterization for
645      !!  ==  transfer coefficients  ==!   Cd, Ch, Ce at T-point, and more...
646      SELECT CASE( nblk )
647
648      CASE( np_NCAR      )
649         CALL turb_ncar    ( rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm,                              &
650            &                zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
651
652      CASE( np_COARE_3p0 )
653         CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, &
654            &                zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,   &
655            &                Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) )
656
657      CASE( np_COARE_3p6 )
658         CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, &
659            &                zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,   &
660            &                Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) )
661
662      CASE( np_ECMWF     )
663         CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl,  &
664            &                zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,   &
665            &                Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) )
666
667      CASE DEFAULT
668         CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' )
669
670      END SELECT
671     
672      IF( iom_use('Cd_oce') )   CALL iom_put("Cd_oce",   zcd_oce * tmask(:,:,1))
673      IF( iom_use('Ce_oce') )   CALL iom_put("Ce_oce",   zce_oce * tmask(:,:,1))
674      IF( iom_use('Ch_oce') )   CALL iom_put("Ch_oce",   zch_oce * tmask(:,:,1))
675      !! LB: mainly here for debugging purpose:
676      IF( iom_use('theta_zt') ) CALL iom_put("theta_zt", (ztpot-rt0) * tmask(:,:,1)) ! potential temperature at z=zt
677      IF( iom_use('q_zt') )     CALL iom_put("q_zt",     zqair       * tmask(:,:,1)) ! specific humidity       "
678      IF( iom_use('theta_zu') ) CALL iom_put("theta_zu", (t_zu -rt0) * tmask(:,:,1)) ! potential temperature at z=zu
679      IF( iom_use('q_zu') )     CALL iom_put("q_zu",     q_zu        * tmask(:,:,1)) ! specific humidity       "
680      IF( iom_use('ssq') )      CALL iom_put("ssq",      pssq        * tmask(:,:,1)) ! saturation specific humidity at z=0
681      IF( iom_use('wspd_blk') ) CALL iom_put("wspd_blk", zU_zu       * tmask(:,:,1)) ! bulk wind speed at z=zu
682     
683      IF( ln_skin_cs .OR. ln_skin_wl ) THEN
684         !! ptsk and pssq have been updated!!!
685         !!
686         !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of ptsk and pssq:
687         WHERE ( fr_i(:,:) > 0.001_wp )
688            ! sea-ice present, we forget about the update, using what we backed up before call to turb_*()
689            ptsk(:,:) = zztmp1(:,:)
690            pssq(:,:) = zztmp2(:,:)
691         END WHERE
692      END IF
693
694      !  Turbulent fluxes over ocean  => BULK_FORMULA @ sbcblk_phy.F90
695      ! -------------------------------------------------------------
696
697      IF( ln_abl ) THEN         !==  ABL formulation  ==!   multiplication by rho_air and turbulent fluxes computation done in ablstp
698         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
699            zztmp = zU_zu(ji,jj)
700            wndm(ji,jj)   = zztmp                   ! Store zU_zu in wndm to compute ustar2 in ablmod
701            pcd_du(ji,jj) = zztmp * zcd_oce(ji,jj)
702            psen(ji,jj)   = zztmp * zch_oce(ji,jj)
703            pevp(ji,jj)   = zztmp * zce_oce(ji,jj)
704            rhoa(ji,jj)   = rho_air( ptair(ji,jj), phumi(ji,jj), pslp(ji,jj) )
705         END_2D
706      ELSE                      !==  BLK formulation  ==!   turbulent fluxes computation
707         CALL BULK_FORMULA( rn_zu, ptsk(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), &
708            &               zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:),          &
709            &               wndm(:,:), zU_zu(:,:), pslp(:,:),                  &
710            &               taum(:,:), psen(:,:), zqla(:,:),                   &
711            &               pEvap=pevp(:,:), prhoa=rhoa(:,:), pfact_evap=rn_efac )
712
713         zqla(:,:) = zqla(:,:) * tmask(:,:,1)
714         psen(:,:) = psen(:,:) * tmask(:,:,1)
715         taum(:,:) = taum(:,:) * tmask(:,:,1)
716         pevp(:,:) = pevp(:,:) * tmask(:,:,1)
717
718         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
719            IF( wndm(ji,jj) > 0._wp ) THEN
720               zztmp = taum(ji,jj) / wndm(ji,jj)
721#if defined key_cyclone
722               ztau_i(ji,jj) = zztmp * zwnd_i(ji,jj)
723               ztau_j(ji,jj) = zztmp * zwnd_j(ji,jj)
724#else
725               ztau_i(ji,jj) = zztmp * pwndi(ji,jj)
726               ztau_j(ji,jj) = zztmp * pwndj(ji,jj)
727#endif
728            ELSE
729               ztau_i(ji,jj) = 0._wp
730               ztau_j(ji,jj) = 0._wp                 
731            ENDIF
732         END_2D
733
734         IF( ln_crt_fbk ) THEN   ! aply eq. 10 and 11 of Renault et al. 2020 (doi: 10.1029/2019MS001715)
735            zstmax = MIN( rn_stau_a * 3._wp + rn_stau_b, 0._wp )   ! set the max value of Stau corresponding to a wind of 3 m/s (<0)
736            DO_2D( 0, 1, 0, 1 )   ! end at jpj and jpi, as ztau_j(ji,jj+1) ztau_i(ji+1,jj) used in the next loop
737               zstau = MIN( rn_stau_a * wndm(ji,jj) + rn_stau_b, zstmax )   ! stau (<0) must be smaller than zstmax
738               ztau_i(ji,jj) = ztau_i(ji,jj) + zstau * ( 0.5_wp * ( pu(ji-1,jj  ) + pu(ji,jj) ) - puatm(ji,jj) )
739               ztau_j(ji,jj) = ztau_j(ji,jj) + zstau * ( 0.5_wp * ( pv(ji  ,jj-1) + pv(ji,jj) ) - pvatm(ji,jj) )
740               taum(ji,jj) = SQRT( ztau_i(ji,jj) * ztau_i(ji,jj) + ztau_j(ji,jj) * ztau_j(ji,jj) )
741            END_2D
742         ENDIF
743
744         ! ... utau, vtau at U- and V_points, resp.
745         !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
746         !     Note that coastal wind stress is not used in the code... so this extra care has no effect
747         DO_2D( 0, 0, 0, 0 )              ! start loop at 2, in case ln_crt_fbk = T
748            utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( ztau_i(ji,jj) + ztau_i(ji+1,jj  ) ) &
749               &              * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))
750            vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( ztau_j(ji,jj) + ztau_j(ji  ,jj+1) ) &
751               &              * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1))
752         END_2D
753
754         IF( ln_crt_fbk ) THEN
755            CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1., taum, 'T', -1. )
756         ELSE
757            CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. )
758         ENDIF
759
760         CALL iom_put( "taum_oce", taum )   ! output wind stress module
761
762         IF(sn_cfctl%l_prtctl) THEN
763            CALL prt_ctl( tab2d_1=wndm  , clinfo1=' blk_oce_1: wndm   : ')
764            CALL prt_ctl( tab2d_1=utau  , clinfo1=' blk_oce_1: utau   : ', mask1=umask,   &
765               &          tab2d_2=vtau  , clinfo2='            vtau   : ', mask2=vmask )
766         ENDIF
767         !
768      ENDIF !IF( ln_abl )
769     
770      ptsk(:,:) = ( ptsk(:,:) - rt0 ) * tmask(:,:,1)  ! Back to Celsius
771           
772      IF( ln_skin_cs .OR. ln_skin_wl ) THEN
773         CALL iom_put( "t_skin" ,  ptsk        )  ! T_skin in Celsius
774         CALL iom_put( "dt_skin" , ptsk - pst  )  ! T_skin - SST temperature difference...
775      ENDIF
776
777      IF(sn_cfctl%l_prtctl) THEN
778         CALL prt_ctl( tab2d_1=pevp  , clinfo1=' blk_oce_1: pevp   : ' )
779         CALL prt_ctl( tab2d_1=psen  , clinfo1=' blk_oce_1: psen   : ' )
780         CALL prt_ctl( tab2d_1=pssq  , clinfo1=' blk_oce_1: pssq   : ' )
781      ENDIF
782      !
783   END SUBROUTINE blk_oce_1
784
785
786   SUBROUTINE blk_oce_2( ptair, pqsr, pqlw, pprec,   &   ! <<= in
787      &                  psnow, ptsk, psen, pevp     )   ! <<= in
788      !!---------------------------------------------------------------------
789      !!                     ***  ROUTINE blk_oce_2  ***
790      !!
791      !! ** Purpose :   finalize the momentum, heat and freshwater fluxes computation
792      !!                at the ocean surface at each time step knowing Cd, Ch, Ce and
793      !!                atmospheric variables (from ABL or external data)
794      !!
795      !! ** Outputs : - utau    : i-component of the stress at U-point  (N/m2)
796      !!              - vtau    : j-component of the stress at V-point  (N/m2)
797      !!              - taum    : Wind stress module at T-point         (N/m2)
798      !!              - wndm    : Wind speed module at T-point          (m/s)
799      !!              - qsr     : Solar heat flux over the ocean        (W/m2)
800      !!              - qns     : Non Solar heat flux over the ocean    (W/m2)
801      !!              - emp     : evaporation minus precipitation       (kg/m2/s)
802      !!---------------------------------------------------------------------
803      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptair
804      REAL(wp), INTENT(in), DIMENSION(:,:) ::   pqsr
805      REAL(wp), INTENT(in), DIMENSION(:,:) ::   pqlw
806      REAL(wp), INTENT(in), DIMENSION(:,:) ::   pprec
807      REAL(wp), INTENT(in), DIMENSION(:,:) ::   psnow
808      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptsk   ! SKIN surface temperature   [Celsius]
809      REAL(wp), INTENT(in), DIMENSION(:,:) ::   psen
810      REAL(wp), INTENT(in), DIMENSION(:,:) ::   pevp
811      !
812      INTEGER  ::   ji, jj               ! dummy loop indices
813      REAL(wp) ::   zztmp,zz1,zz2,zz3    ! local variable
814      REAL(wp), DIMENSION(jpi,jpj) ::   ztskk             ! skin temp. in Kelvin
815      REAL(wp), DIMENSION(jpi,jpj) ::   zqlw              ! long wave and sensible heat fluxes     
816      REAL(wp), DIMENSION(jpi,jpj) ::   zqla              ! latent heat fluxes and evaporation
817      !!---------------------------------------------------------------------
818      !
819      ! local scalars ( place there for vector optimisation purposes)
820
821
822      ztskk(:,:) = ptsk(:,:) + rt0  ! => ptsk in Kelvin rather than Celsius
823     
824      ! ----------------------------------------------------------------------------- !
825      !     III    Net longwave radiative FLUX                                        !
826      ! ----------------------------------------------------------------------------- !
827
828      !! LB: now moved after Turbulent fluxes because must use the skin temperature rather that the SST
829      !! (ztskk is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.)
830      zqlw(:,:) = emiss_w * ( pqlw(:,:) - stefan*ztskk(:,:)*ztskk(:,:)*ztskk(:,:)*ztskk(:,:) ) * tmask(:,:,1)   ! Net radiative longwave flux
831
832      !  Latent flux over ocean
833      ! -----------------------
834
835      ! use scalar version of L_vap() for AGRIF compatibility
836      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
837         zqla(ji,jj) = - L_vap( ztskk(ji,jj) ) * pevp(ji,jj)    ! Latent Heat flux !!GS: possibility to add a global qla to avoid recomputation after abl update
838      END_2D
839
840      IF(sn_cfctl%l_prtctl) THEN
841         CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce_2: zqla   : ' )
842         CALL prt_ctl( tab2d_1=zqlw  , clinfo1=' blk_oce_2: zqlw   : ', tab2d_2=qsr, clinfo2=' qsr : ' )
843
844      ENDIF
845
846      ! ----------------------------------------------------------------------------- !
847      !     IV    Total FLUXES                                                       !
848      ! ----------------------------------------------------------------------------- !
849      !
850      emp (:,:) = (  pevp(:,:)                                       &   ! mass flux (evap. - precip.)
851         &         - pprec(:,:) * rn_pfac  ) * tmask(:,:,1)
852      !
853      qns(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:)                   &   ! Downward Non Solar
854         &     - psnow(:,:) * rn_pfac * rLfus                        &   ! remove latent melting heat for solid precip
855         &     - pevp(:,:) * ptsk(:,:) * rcp                         &   ! remove evap heat content at SST
856         &     + ( pprec(:,:) - psnow(:,:) ) * rn_pfac               &   ! add liquid precip heat content at Tair
857         &     * ( ptair(:,:) - rt0 ) * rcp                          &
858         &     + psnow(:,:) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow)
859         &     * ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi
860      qns(:,:) = qns(:,:) * tmask(:,:,1)
861      !
862#if defined key_si3
863      qns_oce(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:)                             ! non solar without emp (only needed by SI3)
864      qsr_oce(:,:) = qsr(:,:)
865#endif
866      !
867      CALL iom_put( "rho_air"  , rhoa*tmask(:,:,1) )       ! output air density [kg/m^3]
868      CALL iom_put( "evap_oce" , pevp )                    ! evaporation
869      CALL iom_put( "qlw_oce"  , zqlw )                    ! output downward longwave heat over the ocean
870      CALL iom_put( "qsb_oce"  , psen )                    ! output downward sensible heat over the ocean
871      CALL iom_put( "qla_oce"  , zqla )                    ! output downward latent   heat over the ocean
872      tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1)   ! output total precipitation [kg/m2/s]
873      sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1)   ! output solid precipitation [kg/m2/s]
874      CALL iom_put( 'snowpre', sprecip )                   ! Snow
875      CALL iom_put( 'precip' , tprecip )                   ! Total precipitation
876      !
877      IF ( nn_ice == 0 ) THEN
878         CALL iom_put( "qemp_oce" , qns-zqlw-psen-zqla )   ! output downward heat content of E-P over the ocean
879         CALL iom_put( "qns_oce"  ,   qns  )               ! output downward non solar heat over the ocean
880         CALL iom_put( "qsr_oce"  ,   qsr  )               ! output downward solar heat over the ocean
881         CALL iom_put( "qt_oce"   ,   qns+qsr )            ! output total downward heat over the ocean
882      ENDIF
883      !
884      IF(sn_cfctl%l_prtctl) THEN
885         CALL prt_ctl(tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw  : ')
886         CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce_2: zqla  : ', tab2d_2=qsr  , clinfo2=' qsr   : ')
887         CALL prt_ctl(tab2d_1=emp  , clinfo1=' blk_oce_2: emp   : ')
888      ENDIF
889      !
890   END SUBROUTINE blk_oce_2
891
892
893   SUBROUTINE iom_blk_wgt_init( cdname, ld_tmppatch ) 
894#if defined key_iomput
895      use xios
896#endif
897      !!----------------------------------------------------------------------
898      !!                     ***  ROUTINE   ***
899      !!
900      !! ** Purpose :   initialize context for reading weights for surface
901      !!                forcing
902      !!
903      !!----------------------------------------------------------------------
904      CHARACTER(len=*),           INTENT(in)  :: cdname
905      LOGICAL         , OPTIONAL, INTENT(in)  :: ld_tmppatch
906#if defined key_iomput
907      !
908      INTEGER, PARAMETER :: lcname = 26
909      CHARACTER(len=lc) :: clname
910      INTEGER           :: ji, jkmin
911      LOGICAL ::   ll_tmppatch = .TRUE.    !: seb: patch before we remove periodicity
912      INTEGER ::   nldi_save, nlei_save    !:      and close boundaries in output files
913      INTEGER ::   nldj_save, nlej_save    !:
914      LOGICAL ::   ll_global = .FALSE.     !: do we have variable on model grid
915      CHARACTER(len=lc), DIMENSION( jpts) :: cg_name(jpts)
916      CHARACTER(len=1), DIMENSION(lcname) :: cname
917      CHARACTER(len=lc)               :: cfname ! file name without .nc
918      TYPE(xios_duration)               :: dtime    = xios_duration(0, 0, 0, 0, 0, 0), &
919                                           outp_frq = xios_duration(0, 0, 0, 0, 0, 0)
920      TYPE(xios_domaingroup)            :: domaingroup_hdl
921      TYPE(xios_domain)                 :: domain_hdl
922      TYPE(xios_axisgroup)              :: axisgroup_hdl
923      TYPE(xios_axis)                   :: axis_hdl
924      TYPE(xios_scalar)                 :: scalar_hdl
925      TYPE(xios_scalargroup)            :: scalargroup_hdl
926      TYPE(xios_file)                   :: file_hdl
927      TYPE(xios_filegroup)              :: filegroup_hdl
928      TYPE(xios_field)                  :: field_hdl
929      INTEGER                           :: jf, ni, nj, ipos, jfld
930      INTEGER, DIMENSION(3)             :: ndims         ! size of each dimension. this is surface
931                                                         ! forcing, not more than 3: (x, y, t)
932      INTEGER                           :: ndim          ! number if dimensions
933      INTEGER                           :: inum          ! file id
934      INTEGER                           :: ivid          ! varable id in netcdf file
935      CHARACTER (len=5)                 :: aname         !name of the field in wght file
936      INTEGER                           :: numwgt        !
937      LOGICAL                           :: lfound
938      INTEGER                           :: jg, id, jn
939     
940
941      cname(:)= (/'a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', &
942                  'k', 'l', 'm', 'n', 'o', 'p', 'q', 'r', 's', 't', &
943                  'u', 'v', 'w', 'x', 'y', 'z'/)
944      !!----------------------------------------------------------------------
945      !
946      ! seb: patch before we remove periodicity and close boundaries in output files
947      IF ( ll_tmppatch ) THEN
948         nldi_save = nldi   ;   nlei_save = nlei
949         nldj_save = nldj   ;   nlej_save = nlej
950         IF( nimpp           ==      1 ) nldi = 1
951         IF( nimpp + jpi - 1 == jpiglo ) nlei = jpi
952         IF( njmpp           ==      1 ) nldj = 1
953         IF( njmpp + jpj - 1 == jpjglo ) nlej = jpj
954      ENDIF
955      !
956      jfld = SIZE(sf)
957
958      IF(jfld > lcname) THEN
959      STOP 'lcname number is too small'
960      ENDIF
961
962      clname = cdname
963      IF( TRIM(Agrif_CFixed()) /= '0' )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(cdname)
964      CALL xios_context_initialize(TRIM(clname), mpi_comm_oce)
965      CALL iom_swap( cdname )
966      ! Calendar type is now defined in xml file
967      SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
968      CASE ( 1)   ; CALL xios_define_calendar( TYPE = "Gregorian", time_origin = xios_date(1948, 01, 01, 00, 00, 00), &
969          &                                    start_date = xios_date(nyear, nmonth, nday,00,00,00) )
970      CASE ( 0)   ; CALL xios_define_calendar( TYPE = "NoLeap"   , time_origin = xios_date(1948, 01, 01, 00, 00, 00), &
971          &                                    start_date = xios_date(nyear, nmonth, nday,00,00,00) )
972      CASE (30)   ; CALL xios_define_calendar( TYPE = "D360"     , time_origin = xios_date(1948, 01, 01, 00, 00, 00), &
973          &                                    start_date = xios_date(nyear, nmonth, nday, 00, 00, 00) )
974      END SELECT
975
976
977      dtime%month = 1
978      CALL xios_set_timestep( dtime )
979
980      CALL xios_get_handle("domain_definition",domaingroup_hdl)
981      CALL xios_add_child(domaingroup_hdl, domain_hdl, "grid_global")
982      ni = nlei-nldi+1
983      nj = nlej-nldj+1
984      CALL xios_set_domain_attr("grid_global", ni_glo=jpiglo, nj_glo=jpjglo, ibegin=nimpp+nldi-2, jbegin=njmpp+nldj-2, ni=ni, nj=nj)
985      CALL xios_set_domain_attr("grid_global", data_dim=2, data_ibegin = 1-nldi, data_ni = jpi, data_jbegin = 1-nldj, data_nj = jpj)
986      CALL xios_set_domain_attr("grid_global", type='curvilinear')
987
988      CALL xios_get_handle("file_definition", filegroup_hdl )
989
990      DO jf = 1, jfld
991!do we have weights file defined already?
992         lfound = .FALSE.
993         DO jg = 1, jf-1
994            IF(TRIM(sf(jf)%wgtname) == TRIM(sf(jg)%wgtname)) THEN
995               IF(lwp) write(numout, *) TRIM(sf(jf)%wgtname),' already defined, skipping '
996               lfound = .TRUE.
997               exit
998            ENDIF
999         ENDDO
1000         IF(lfound) cycle
1001!add next weights file
1002         ipos = index(sf(jf)%wgtname,'.nc')
1003         cfname(1:lc) = " "
1004         IF(ipos > 0) THEN
1005          cfname(1:ipos-1) = sf(jf)%wgtname(1:ipos-1)
1006         ELSE
1007          cfname(1:lc) = sf(jf)%wgtname(1:lc)
1008         ENDIF
1009
1010         CALL xios_add_child(filegroup_hdl, file_hdl, cname(jf))
1011         CALL xios_set_file_attr( cname(jf), name=TRIM(cfname),           &
1012                type="one_file", par_access="collective", enabled=.TRUE., &
1013                output_freq=xios_timestep, mode="read")
1014
1015         CALL iom_open ( sf(jf)%wgtname, inum )
1016         id = iom_varid(inum, 'src05', ldstop=.FALSE.)
1017         IF(id <= 0) THEN
1018             numwgt = 4
1019         ELSE
1020             numwgt = 16
1021         ENDIF
1022         CALL iom_close ( inum )
1023!define variables in weights file. Because all files have the same names src01,
1024!dst01, ... add filename to the name to uniquely identify variables
1025         aname = ' '
1026         DO jn = 1, numwgt
1027            WRITE(aname,'(a3,i2.2)') 'src',jn
1028            CALL xios_add_child(file_hdl, field_hdl, TRIM(aname)//TRIM(sf(jf)%wgtname))
1029            CALL xios_set_attr (field_hdl, enabled = .TRUE., &
1030                          name = TRIM(aname), domain_ref="grid_global", &
1031                          operation = "instant")
1032            WRITE(aname,'(a3,i2.2)') 'dst',jn
1033            CALL xios_add_child(file_hdl, field_hdl, TRIM(aname)//TRIM(sf(jf)%wgtname))
1034            CALL xios_set_attr (field_hdl, enabled = .TRUE., &
1035                          name = TRIM(aname), domain_ref="grid_global", &
1036                          operation = "instant")
1037            WRITE(aname,'(a3,i2.2)') 'wgt',jn
1038            CALL xios_add_child(file_hdl, field_hdl, TRIM(aname)//TRIM(sf(jf)%wgtname))
1039            CALL xios_set_attr (field_hdl, enabled = .TRUE., &
1040                          name = TRIM(aname), domain_ref="grid_global", &
1041                          operation = "instant")
1042         ENDDO 
1043      ENDDO   
1044
1045
1046      CALL xios_close_context_definition()
1047      CALL xios_update_calendar( 0 )
1048
1049      IF ( ll_tmppatch ) THEN
1050         nldi = nldi_save   ;   nlei = nlei_save
1051         nldj = nldj_save   ;   nlej = nlej_save
1052      ENDIF
1053#endif
1054      !
1055   END SUBROUTINE iom_blk_wgt_init
1056
1057
1058#if defined key_si3
1059   !!----------------------------------------------------------------------
1060   !!   'key_si3'                                       SI3 sea-ice model
1061   !!----------------------------------------------------------------------
1062   !!   blk_ice_1   : provide the air-ice stress
1063   !!   blk_ice_2   : provide the heat and mass fluxes at air-ice interface
1064   !!   blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
1065   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
1066   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag
1067   !!----------------------------------------------------------------------
1068
1069   SUBROUTINE blk_ice_1( pwndi, pwndj, ptair, phumi, pslp , puice, pvice, ptsui,  &   ! inputs
1070      &                  putaui, pvtaui, pseni, pevpi, pssqi, pcd_dui             )   ! optional outputs
1071      !!---------------------------------------------------------------------
1072      !!                     ***  ROUTINE blk_ice_1  ***
1073      !!
1074      !! ** Purpose :   provide the surface boundary condition over sea-ice
1075      !!
1076      !! ** Method  :   compute momentum using bulk formulation
1077      !!                formulea, ice variables and read atmospheric fields.
1078      !!                NB: ice drag coefficient is assumed to be a constant
1079      !!---------------------------------------------------------------------
1080      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pslp    ! sea-level pressure [Pa]
1081      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pwndi   ! atmospheric wind at T-point [m/s]
1082      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pwndj   ! atmospheric wind at T-point [m/s]
1083      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   ptair   ! atmospheric wind at T-point [m/s]
1084      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   phumi   ! atmospheric wind at T-point [m/s]
1085      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   puice   ! sea-ice velocity on I or C grid [m/s]
1086      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pvice   ! "
1087      REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   ptsui   ! sea-ice surface temperature [K]
1088      REAL(wp) , INTENT(  out), DIMENSION(:,:  ), OPTIONAL ::   putaui  ! if ln_blk
1089      REAL(wp) , INTENT(  out), DIMENSION(:,:  ), OPTIONAL ::   pvtaui  ! if ln_blk
1090      REAL(wp) , INTENT(  out), DIMENSION(:,:  ), OPTIONAL ::   pseni   ! if ln_abl
1091      REAL(wp) , INTENT(  out), DIMENSION(:,:  ), OPTIONAL ::   pevpi   ! if ln_abl
1092      REAL(wp) , INTENT(  out), DIMENSION(:,:  ), OPTIONAL ::   pssqi   ! if ln_abl
1093      REAL(wp) , INTENT(  out), DIMENSION(:,:  ), OPTIONAL ::   pcd_dui ! if ln_abl
1094      !
1095      INTEGER  ::   ji, jj    ! dummy loop indices
1096      REAL(wp) ::   zootm_su                      ! sea-ice surface mean temperature
1097      REAL(wp) ::   zztmp1, zztmp2                ! temporary arrays
1098      REAL(wp), DIMENSION(jpi,jpj) ::   zcd_dui   ! transfer coefficient for momentum      (tau)
1099      !!---------------------------------------------------------------------
1100      !
1101
1102      ! ------------------------------------------------------------ !
1103      !    Wind module relative to the moving ice ( U10m - U_ice )   !
1104      ! ------------------------------------------------------------ !
1105      ! C-grid ice dynamics :   U & V-points (same as ocean)
1106      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
1107         wndm_ice(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) )
1108      END_2D
1109      !
1110      ! Make ice-atm. drag dependent on ice concentration
1111      IF    ( ln_Cd_L12 ) THEN   ! calculate new drag from Lupkes(2012) equations
1112         CALL Cdn10_Lupkes2012( Cd_ice )
1113         Ch_ice(:,:) = Cd_ice(:,:)       ! momentum and heat transfer coef. are considered identical
1114         Ce_ice(:,:) = Cd_ice(:,:)
1115      ELSEIF( ln_Cd_L15 ) THEN   ! calculate new drag from Lupkes(2015) equations
1116         CALL Cdn10_Lupkes2015( ptsui, pslp, Cd_ice, Ch_ice )
1117         Ce_ice(:,:) = Ch_ice(:,:)       ! sensible and latent heat transfer coef. are considered identical
1118      ENDIF
1119     
1120      IF( iom_use('Cd_ice') ) CALL iom_put("Cd_ice", Cd_ice)
1121      IF( iom_use('Ce_ice') ) CALL iom_put("Ce_ice", Ce_ice)
1122      IF( iom_use('Ch_ice') ) CALL iom_put("Ch_ice", Ch_ice)
1123     
1124      ! local scalars ( place there for vector optimisation purposes)
1125      zcd_dui(:,:) = wndm_ice(:,:) * Cd_ice(:,:)
1126
1127      IF( ln_blk ) THEN
1128         ! ---------------------------------------------------- !
1129         !    Wind stress relative to nonmoving ice ( U10m )    !
1130         ! ---------------------------------------------------- !
1131         ! supress moving ice in wind stress computation as we don't know how to do it properly...
1132         DO_2D( 0, 1, 0, 1 )    ! at T point
1133            putaui(ji,jj) = rhoa(ji,jj) * zcd_dui(ji,jj) * pwndi(ji,jj)
1134            pvtaui(ji,jj) = rhoa(ji,jj) * zcd_dui(ji,jj) * pwndj(ji,jj)
1135         END_2D
1136         !
1137         DO_2D( 0, 0, 0, 0 )    ! U & V-points (same as ocean).
1138            ! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and  rheology
1139            zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj  ,1) )
1140            zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji  ,jj+1,1) )
1141            putaui(ji,jj) = zztmp1 * ( putaui(ji,jj) + putaui(ji+1,jj  ) )
1142            pvtaui(ji,jj) = zztmp2 * ( pvtaui(ji,jj) + pvtaui(ji  ,jj+1) )
1143         END_2D
1144         CALL lbc_lnk_multi( 'sbcblk', putaui, 'U', -1._wp, pvtaui, 'V', -1._wp )
1145         !
1146         IF(sn_cfctl%l_prtctl)  CALL prt_ctl( tab2d_1=putaui  , clinfo1=' blk_ice: putaui : '   &
1147            &                               , tab2d_2=pvtaui  , clinfo2='          pvtaui : ' )
1148      ELSE ! ln_abl
1149         zztmp1 = 11637800.0_wp
1150         zztmp2 =    -5897.8_wp
1151         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
1152            pcd_dui(ji,jj) = zcd_dui (ji,jj)
1153            pseni  (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj)
1154            pevpi  (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj)
1155            zootm_su       = zztmp2 / ptsui(ji,jj)   ! ptsui is in K (it can't be zero ??)
1156            pssqi  (ji,jj) = zztmp1 * EXP( zootm_su ) / rhoa(ji,jj)
1157         END_2D
1158      ENDIF
1159      !
1160      IF(sn_cfctl%l_prtctl)  CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice: wndm_ice : ')
1161      !
1162   END SUBROUTINE blk_ice_1
1163
1164
1165   SUBROUTINE blk_ice_2( ptsu, phs, phi, palb, ptair, phumi, pslp, pqlw, pprec, psnow  )
1166      !!---------------------------------------------------------------------
1167      !!                     ***  ROUTINE blk_ice_2  ***
1168      !!
1169      !! ** Purpose :   provide the heat and mass fluxes at air-ice interface
1170      !!
1171      !! ** Method  :   compute heat and freshwater exchanged
1172      !!                between atmosphere and sea-ice using bulk formulation
1173      !!                formulea, ice variables and read atmmospheric fields.
1174      !!
1175      !! caution : the net upward water flux has with mm/day unit
1176      !!---------------------------------------------------------------------
1177      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu   ! sea ice surface temperature [K]
1178      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness
1179      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness
1180      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb   ! ice albedo (all skies)
1181      REAL(wp), DIMENSION(:,:  ), INTENT(in)  ::   ptair
1182      REAL(wp), DIMENSION(:,:  ), INTENT(in)  ::   phumi
1183      REAL(wp), DIMENSION(:,:  ), INTENT(in)  ::   pslp
1184      REAL(wp), DIMENSION(:,:  ), INTENT(in)  ::   pqlw
1185      REAL(wp), DIMENSION(:,:  ), INTENT(in)  ::   pprec
1186      REAL(wp), DIMENSION(:,:  ), INTENT(in)  ::   psnow
1187      !!
1188      INTEGER  ::   ji, jj, jl               ! dummy loop indices
1189      REAL(wp) ::   zst3                     ! local variable
1190      REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      -
1191      REAL(wp) ::   zztmp, zztmp2, z1_rLsub  !   -      -
1192      REAL(wp) ::   zfr1, zfr2               ! local variables
1193      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_st         ! inverse of surface temperature
1194      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qlw         ! long wave heat flux over ice
1195      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qsb         ! sensible  heat flux over ice
1196      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_dqlw        ! long wave heat sensitivity over ice
1197      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_dqsb        ! sensible  heat sensitivity over ice
1198      REAL(wp), DIMENSION(jpi,jpj)     ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (SI3)
1199      REAL(wp), DIMENSION(jpi,jpj)     ::   zqair         ! specific humidity of air at z=rn_zqt [kg/kg] !LB
1200      REAL(wp), DIMENSION(jpi,jpj)     ::   ztmp, ztmp2
1201      !!---------------------------------------------------------------------
1202      !
1203      zcoef_dqlw = 4._wp * 0.95_wp * stefan             ! local scalars
1204      zcoef_dqla = -rLsub * 11637800._wp * (-5897.8_wp) !LB: BAD!
1205      !
1206      SELECT CASE( nhumi )
1207      CASE( np_humi_sph )
1208         zqair(:,:) =  phumi(:,:)      ! what we read in file is already a spec. humidity!
1209      CASE( np_humi_dpt )
1210         zqair(:,:) = q_sat( phumi(:,:), pslp )
1211      CASE( np_humi_rlh )
1212         zqair(:,:) = q_air_rh( 0.01_wp*phumi(:,:), ptair(:,:), pslp(:,:) ) !LB: 0.01 => RH is % percent in file
1213      END SELECT
1214      !
1215      zztmp = 1. / ( 1. - albo )
1216      WHERE( ptsu(:,:,:) /= 0._wp )
1217         z1_st(:,:,:) = 1._wp / ptsu(:,:,:)
1218      ELSEWHERE
1219         z1_st(:,:,:) = 0._wp
1220      END WHERE
1221      !                                     ! ========================== !
1222      DO jl = 1, jpl                        !  Loop over ice categories  !
1223         !                                  ! ========================== !
1224         DO jj = 1 , jpj
1225            DO ji = 1, jpi
1226               ! ----------------------------!
1227               !      I   Radiative FLUXES   !
1228               ! ----------------------------!
1229               zst3 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl) * ptsu(ji,jj,jl)
1230               ! Short Wave (sw)
1231               qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
1232               ! Long  Wave (lw)
1233               z_qlw(ji,jj,jl) = 0.95 * ( pqlw(ji,jj) - stefan * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
1234               ! lw sensitivity
1235               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
1236
1237               ! ----------------------------!
1238               !     II    Turbulent FLUXES  !
1239               ! ----------------------------!
1240
1241               ! ... turbulent heat fluxes with Ch_ice recalculated in blk_ice_1
1242               ! Sensible Heat
1243               z_qsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_ice(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - ptair(ji,jj))
1244               ! Latent Heat
1245               zztmp2 = EXP( -5897.8 * z1_st(ji,jj,jl) )
1246               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa(ji,jj) * rLsub  * Ce_ice(ji,jj) * wndm_ice(ji,jj) *  &
1247                  &                ( 11637800. * zztmp2 / rhoa(ji,jj) - zqair(ji,jj) ) )
1248               ! Latent heat sensitivity for ice (Dqla/Dt)
1249               IF( qla_ice(ji,jj,jl) > 0._wp ) THEN
1250                  dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ce_ice(ji,jj) * wndm_ice(ji,jj) *  &
1251                     &                 z1_st(ji,jj,jl) * z1_st(ji,jj,jl) * zztmp2
1252               ELSE
1253                  dqla_ice(ji,jj,jl) = 0._wp
1254               ENDIF
1255
1256               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice)
1257               z_dqsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_ice(ji,jj) * wndm_ice(ji,jj)
1258
1259               ! ----------------------------!
1260               !     III    Total FLUXES     !
1261               ! ----------------------------!
1262               ! Downward Non Solar flux
1263               qns_ice (ji,jj,jl) =     z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl)
1264               ! Total non solar heat flux sensitivity for ice
1265               dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) )
1266            END DO
1267            !
1268         END DO
1269         !
1270      END DO
1271      !
1272      tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1)  ! total precipitation [kg/m2/s]
1273      sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1)  ! solid precipitation [kg/m2/s]
1274      CALL iom_put( 'snowpre', sprecip )                  ! Snow precipitation
1275      CALL iom_put( 'precip' , tprecip )                  ! Total precipitation
1276
1277      ! --- evaporation --- !
1278      z1_rLsub = 1._wp / rLsub
1279      evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_rLsub    ! sublimation
1280      devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_rLsub    ! d(sublimation)/dT
1281      zevap    (:,:)   = emp(:,:) + tprecip(:,:)   ! evaporation over ocean  !LB: removed rn_efac here, correct???
1282
1283      ! --- evaporation minus precipitation --- !
1284      zsnw(:,:) = 0._wp
1285      CALL ice_thd_snwblow( (1.-at_i_b(:,:)), zsnw )  ! snow distribution over ice after wind blowing
1286      emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
1287      emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
1288      emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:)
1289
1290      ! --- heat flux associated with emp --- !
1291      qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp                  & ! evap at sst
1292         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( ptair(:,:) - rt0 ) * rcp               & ! liquid precip at Tair
1293         &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip at min(Tair,Tsnow)
1294         &              ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
1295      qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only)
1296         &              ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
1297
1298      ! --- total solar and non solar fluxes --- !
1299      qns_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 )  &
1300         &           + qemp_ice(:,:) + qemp_oce(:,:)
1301      qsr_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
1302
1303      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
1304      qprec_ice(:,:) = rhos * ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
1305
1306      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) ---
1307      DO jl = 1, jpl
1308         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * rcpi * tmask(:,:,1) )
1309         !                         ! But we do not have Tice => consider it at 0degC => evap=0
1310      END DO
1311
1312      ! --- shortwave radiation transmitted below the surface (W/m2, see Grenfell Maykut 77) --- !
1313      zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            ! transmission when hi>10cm
1314      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            ! zfr2 such that zfr1 + zfr2 to equal 1
1315      !
1316      WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm
1317         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) )
1318      ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (zfr1) when hi>10cm
1319         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1
1320      ELSEWHERE                                                         ! zero when hs>0
1321         qtr_ice_top(:,:,:) = 0._wp
1322      END WHERE
1323      !
1324
1325      IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') ) THEN
1326         ztmp(:,:) = zevap(:,:) * ( 1._wp - at_i_b(:,:) )
1327         IF( iom_use('evap_ao_cea'  ) )  CALL iom_put( 'evap_ao_cea'  , ztmp(:,:) * tmask(:,:,1) )   ! ice-free oce evap (cell average)
1328         IF( iom_use('hflx_evap_cea') )  CALL iom_put( 'hflx_evap_cea', ztmp(:,:) * sst_m(:,:) * rcp * tmask(:,:,1) )   ! heat flux from evap (cell average)
1329      ENDIF
1330      IF( iom_use('hflx_rain_cea') ) THEN
1331         ztmp(:,:) = rcp * ( SUM( (ptsu-rt0) * a_i_b, dim=3 ) + sst_m(:,:) * ( 1._wp - at_i_b(:,:) ) )
1332         IF( iom_use('hflx_rain_cea') )  CALL iom_put( 'hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * ztmp(:,:) )   ! heat flux from rain (cell average)
1333      ENDIF
1334      IF( iom_use('hflx_snow_cea') .OR. iom_use('hflx_snow_ao_cea') .OR. iom_use('hflx_snow_ai_cea')  )  THEN
1335         WHERE( SUM( a_i_b, dim=3 ) > 1.e-10 )
1336            ztmp(:,:) = rcpi * SUM( (ptsu-rt0) * a_i_b, dim=3 ) / SUM( a_i_b, dim=3 )
1337         ELSEWHERE
1338            ztmp(:,:) = rcp * sst_m(:,:)
1339         ENDWHERE
1340         ztmp2(:,:) = sprecip(:,:) * ( ztmp(:,:) - rLfus )
1341         IF( iom_use('hflx_snow_cea')    ) CALL iom_put('hflx_snow_cea'   , ztmp2(:,:) ) ! heat flux from snow (cell average)
1342         IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', ztmp2(:,:) * ( 1._wp - zsnw(:,:) ) ) ! heat flux from snow (over ocean)
1343         IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', ztmp2(:,:) *           zsnw(:,:)   ) ! heat flux from snow (over ice)
1344      ENDIF
1345      !
1346      IF(sn_cfctl%l_prtctl) THEN
1347         CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb    : ', kdim=jpl)
1348         CALL prt_ctl(tab3d_1=z_qlw   , clinfo1=' blk_ice: z_qlw    : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl)
1349         CALL prt_ctl(tab3d_1=z_dqsb  , clinfo1=' blk_ice: z_dqsb   : ', tab3d_2=z_dqlw  , clinfo2=' z_dqlw   : ', kdim=jpl)
1350         CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice: dqns_ice : ', tab3d_2=qsr_ice , clinfo2=' qsr_ice  : ', kdim=jpl)
1351         CALL prt_ctl(tab3d_1=ptsu    , clinfo1=' blk_ice: ptsu     : ', tab3d_2=qns_ice , clinfo2=' qns_ice  : ', kdim=jpl)
1352         CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip  : ', tab2d_2=sprecip , clinfo2=' sprecip  : ')
1353      ENDIF
1354      !
1355   END SUBROUTINE blk_ice_2
1356
1357
1358   SUBROUTINE blk_ice_qcn( ld_virtual_itd, ptsu, ptb, phs, phi )
1359      !!---------------------------------------------------------------------
1360      !!                     ***  ROUTINE blk_ice_qcn  ***
1361      !!
1362      !! ** Purpose :   Compute surface temperature and snow/ice conduction flux
1363      !!                to force sea ice / snow thermodynamics
1364      !!                in the case conduction flux is emulated
1365      !!
1366      !! ** Method  :   compute surface energy balance assuming neglecting heat storage
1367      !!                following the 0-layer Semtner (1976) approach
1368      !!
1369      !! ** Outputs : - ptsu    : sea-ice / snow surface temperature (K)
1370      !!              - qcn_ice : surface inner conduction flux (W/m2)
1371      !!
1372      !!---------------------------------------------------------------------
1373      LOGICAL                   , INTENT(in   ) ::   ld_virtual_itd  ! single-category option
1374      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptsu            ! sea ice / snow surface temperature
1375      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   ptb             ! sea ice base temperature
1376      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   phs             ! snow thickness
1377      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   phi             ! sea ice thickness
1378      !
1379      INTEGER , PARAMETER ::   nit = 10                  ! number of iterations
1380      REAL(wp), PARAMETER ::   zepsilon = 0.1_wp         ! characteristic thickness for enhanced conduction
1381      !
1382      INTEGER  ::   ji, jj, jl           ! dummy loop indices
1383      INTEGER  ::   iter                 ! local integer
1384      REAL(wp) ::   zfac, zfac2, zfac3   ! local scalars
1385      REAL(wp) ::   zkeff_h, ztsu, ztsu0 !
1386      REAL(wp) ::   zqc, zqnet           !
1387      REAL(wp) ::   zhe, zqa0            !
1388      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zgfac   ! enhanced conduction factor
1389      !!---------------------------------------------------------------------
1390
1391      ! -------------------------------------!
1392      !      I   Enhanced conduction factor  !
1393      ! -------------------------------------!
1394      ! Emulates the enhancement of conduction by unresolved thin ice (ld_virtual_itd = T)
1395      ! Fichefet and Morales Maqueda, JGR 1997
1396      !
1397      zgfac(:,:,:) = 1._wp
1398
1399      IF( ld_virtual_itd ) THEN
1400         !
1401         zfac  = 1._wp /  ( rn_cnd_s + rcnd_i )
1402         zfac2 = EXP(1._wp) * 0.5_wp * zepsilon
1403         zfac3 = 2._wp / zepsilon
1404         !
1405         DO jl = 1, jpl
1406            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
1407               zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac                            ! Effective thickness
1408               IF( zhe >=  zfac2 )   zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor
1409            END_2D
1410         END DO
1411         !
1412      ENDIF
1413
1414      ! -------------------------------------------------------------!
1415      !      II   Surface temperature and conduction flux            !
1416      ! -------------------------------------------------------------!
1417      !
1418      zfac = rcnd_i * rn_cnd_s
1419      !
1420      DO jl = 1, jpl
1421         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
1422            !
1423            zkeff_h = zfac * zgfac(ji,jj,jl) / &                                    ! Effective conductivity of the snow-ice system divided by thickness
1424               &      ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) )
1425            ztsu    = ptsu(ji,jj,jl)                                                ! Store current iteration temperature
1426            ztsu0   = ptsu(ji,jj,jl)                                                ! Store initial surface temperature
1427            zqa0    = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux
1428            !
1429            DO iter = 1, nit     ! --- Iterative loop
1430               zqc   = zkeff_h * ( ztsu - ptb(ji,jj) )                              ! Conduction heat flux through snow-ice system (>0 downwards)
1431               zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc  ! Surface energy budget
1432               ztsu  = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h )              ! Temperature update
1433            END DO
1434            !
1435            ptsu   (ji,jj,jl) = MIN( rt0, ztsu )
1436            qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) )
1437            qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 )
1438            qml_ice(ji,jj,jl) = ( qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) - qcn_ice(ji,jj,jl) )  &
1439               &   * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) )
1440
1441            ! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- !
1442            hfx_err_dif(ji,jj) = hfx_err_dif(ji,jj) - ( dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) ) * a_i_b(ji,jj,jl)
1443
1444         END_2D
1445         !
1446      END DO
1447      !
1448   END SUBROUTINE blk_ice_qcn
1449
1450
1451   SUBROUTINE Cdn10_Lupkes2012( pcd )
1452      !!----------------------------------------------------------------------
1453      !!                      ***  ROUTINE  Cdn10_Lupkes2012  ***
1454      !!
1455      !! ** Purpose :    Recompute the neutral air-ice drag referenced at 10m
1456      !!                 to make it dependent on edges at leads, melt ponds and flows.
1457      !!                 After some approximations, this can be resumed to a dependency
1458      !!                 on ice concentration.
1459      !!
1460      !! ** Method :     The parameterization is taken from Lupkes et al. (2012) eq.(50)
1461      !!                 with the highest level of approximation: level4, eq.(59)
1462      !!                 The generic drag over a cell partly covered by ice can be re-written as follows:
1463      !!
1464      !!                 Cd = Cdw * (1-A) + Cdi * A + Ce * (1-A)**(nu+1/(10*beta)) * A**mu
1465      !!
1466      !!                    Ce = 2.23e-3       , as suggested by Lupkes (eq. 59)
1467      !!                    nu = mu = beta = 1 , as suggested by Lupkes (eq. 59)
1468      !!                    A is the concentration of ice minus melt ponds (if any)
1469      !!
1470      !!                 This new drag has a parabolic shape (as a function of A) starting at
1471      !!                 Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5
1472      !!                 and going down to Cdi(say 1.4e-3) for A=1
1473      !!
1474      !!                 It is theoretically applicable to all ice conditions (not only MIZ)
1475      !!                 => see Lupkes et al (2013)
1476      !!
1477      !! ** References : Lupkes et al. JGR 2012 (theory)
1478      !!                 Lupkes et al. GRL 2013 (application to GCM)
1479      !!
1480      !!----------------------------------------------------------------------
1481      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pcd
1482      REAL(wp), PARAMETER ::   zCe   = 2.23e-03_wp
1483      REAL(wp), PARAMETER ::   znu   = 1._wp
1484      REAL(wp), PARAMETER ::   zmu   = 1._wp
1485      REAL(wp), PARAMETER ::   zbeta = 1._wp
1486      REAL(wp)            ::   zcoef
1487      !!----------------------------------------------------------------------
1488      zcoef = znu + 1._wp / ( 10._wp * zbeta )
1489
1490      ! generic drag over a cell partly covered by ice
1491      !!Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) +  &                        ! pure ocean drag
1492      !!   &      Cd_ice      *           at_i_b(:,:)   +  &                        ! pure ice drag
1493      !!   &      zCe         * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu   ! change due to sea-ice morphology
1494
1495      ! ice-atm drag
1496      pcd(:,:) = rCd_ice +  &                                                         ! pure ice drag
1497         &      zCe     * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp)  ! change due to sea-ice morphology
1498
1499   END SUBROUTINE Cdn10_Lupkes2012
1500
1501
1502   SUBROUTINE Cdn10_Lupkes2015( ptm_su, pslp, pcd, pch )
1503      !!----------------------------------------------------------------------
1504      !!                      ***  ROUTINE  Cdn10_Lupkes2015  ***
1505      !!
1506      !! ** pUrpose :    Alternative turbulent transfert coefficients formulation
1507      !!                 between sea-ice and atmosphere with distinct momentum
1508      !!                 and heat coefficients depending on sea-ice concentration
1509      !!                 and atmospheric stability (no meltponds effect for now).
1510      !!
1511      !! ** Method :     The parameterization is adapted from Lupkes et al. (2015)
1512      !!                 and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme,
1513      !!                 it considers specific skin and form drags (Andreas et al. 2010)
1514      !!                 to compute neutral transfert coefficients for both heat and
1515      !!                 momemtum fluxes. Atmospheric stability effect on transfert
1516      !!                 coefficient is also taken into account following Louis (1979).
1517      !!
1518      !! ** References : Lupkes et al. JGR 2015 (theory)
1519      !!                 Lupkes et al. ECHAM6 documentation 2015 (implementation)
1520      !!
1521      !!----------------------------------------------------------------------
1522      !
1523      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   ptm_su ! sea-ice surface temperature [K]
1524      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pslp   ! sea-level pressure [Pa]
1525      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pcd    ! momentum transfert coefficient
1526      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pch    ! heat transfert coefficient
1527      REAL(wp), DIMENSION(jpi,jpj)            ::   zst, zqo_sat, zqi_sat
1528      !
1529      ! ECHAM6 constants
1530      REAL(wp), PARAMETER ::   z0_skin_ice  = 0.69e-3_wp  ! Eq. 43 [m]
1531      REAL(wp), PARAMETER ::   z0_form_ice  = 0.57e-3_wp  ! Eq. 42 [m]
1532      REAL(wp), PARAMETER ::   z0_ice       = 1.00e-3_wp  ! Eq. 15 [m]
1533      REAL(wp), PARAMETER ::   zce10        = 2.80e-3_wp  ! Eq. 41
1534      REAL(wp), PARAMETER ::   zbeta        = 1.1_wp      ! Eq. 41
1535      REAL(wp), PARAMETER ::   zc           = 5._wp       ! Eq. 13
1536      REAL(wp), PARAMETER ::   zc2          = zc * zc
1537      REAL(wp), PARAMETER ::   zam          = 2. * zc     ! Eq. 14
1538      REAL(wp), PARAMETER ::   zah          = 3. * zc     ! Eq. 30
1539      REAL(wp), PARAMETER ::   z1_alpha     = 1._wp / 0.2_wp  ! Eq. 51
1540      REAL(wp), PARAMETER ::   z1_alphaf    = z1_alpha    ! Eq. 56
1541      REAL(wp), PARAMETER ::   zbetah       = 1.e-3_wp    ! Eq. 26
1542      REAL(wp), PARAMETER ::   zgamma       = 1.25_wp     ! Eq. 26
1543      REAL(wp), PARAMETER ::   z1_gamma     = 1._wp / zgamma
1544      REAL(wp), PARAMETER ::   r1_3         = 1._wp / 3._wp
1545      !
1546      INTEGER  ::   ji, jj         ! dummy loop indices
1547      REAL(wp) ::   zthetav_os, zthetav_is, zthetav_zu
1548      REAL(wp) ::   zrib_o, zrib_i
1549      REAL(wp) ::   zCdn_skin_ice, zCdn_form_ice, zCdn_ice
1550      REAL(wp) ::   zChn_skin_ice, zChn_form_ice
1551      REAL(wp) ::   z0w, z0i, zfmi, zfmw, zfhi, zfhw
1552      REAL(wp) ::   zCdn_form_tmp
1553      !!----------------------------------------------------------------------
1554
1555      ! Momentum Neutral Transfert Coefficients (should be a constant)
1556      zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2   ! Eq. 40
1557      zCdn_skin_ice = ( vkarmn                                      / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2   ! Eq. 7
1558      zCdn_ice      = zCdn_skin_ice   ! Eq. 7
1559      !zCdn_ice     = 1.89e-3         ! old ECHAM5 value (cf Eq. 32)
1560
1561      ! Heat Neutral Transfert Coefficients
1562      zChn_skin_ice = vkarmn**2 / ( LOG( rn_zu / z0_ice + 1._wp ) * LOG( rn_zu * z1_alpha / z0_skin_ice + 1._wp ) )   ! Eq. 50 + Eq. 52
1563
1564      ! Atmospheric and Surface Variables
1565      zst(:,:)     = sst_m(:,:) + rt0                                        ! convert SST from Celcius to Kelvin
1566      zqo_sat(:,:) = rdct_qsat_salt * q_sat( zst(:,:)   , pslp(:,:) )   ! saturation humidity over ocean [kg/kg]
1567      zqi_sat(:,:) =                  q_sat( ptm_su(:,:), pslp(:,:) )   ! saturation humidity over ice   [kg/kg]
1568      !
1569      DO_2D( 0, 0, 0, 0 )
1570         ! Virtual potential temperature [K]
1571         zthetav_os = zst(ji,jj)    * ( 1._wp + rctv0 * zqo_sat(ji,jj) )   ! over ocean
1572         zthetav_is = ptm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) )   ! ocean ice
1573         zthetav_zu = t_zu (ji,jj)  * ( 1._wp + rctv0 * q_zu(ji,jj)    )   ! at zu
1574
1575         ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead)
1576         zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj)     )**2   ! over ocean
1577         zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2   ! over ice
1578
1579         ! Momentum and Heat Neutral Transfert Coefficients
1580         zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta  ! Eq. 40
1581         zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) )               ! Eq. 53
1582
1583         ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead ?)
1584         z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water
1585         z0i = z0_skin_ice                                             ! over ice
1586         IF( zrib_o <= 0._wp ) THEN
1587            zfmw = 1._wp - zam * zrib_o / ( 1._wp + 3._wp * zc2 * Cdn_oce(ji,jj) * SQRT( -zrib_o * ( rn_zu / z0w + 1._wp ) ) )  ! Eq. 10
1588            zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )**r1_3 / ( Chn_oce(ji,jj) * MAX(0.01, wndm(ji,jj)) )   &     ! Eq. 26
1589               &             )**zgamma )**z1_gamma
1590         ELSE
1591            zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 12
1592            zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 28
1593         ENDIF
1594
1595         IF( zrib_i <= 0._wp ) THEN
1596            zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq.  9
1597            zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq. 25
1598         ELSE
1599            zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 11
1600            zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 27
1601         ENDIF
1602
1603         ! Momentum Transfert Coefficients (Eq. 38)
1604         pcd(ji,jj) = zCdn_skin_ice *   zfmi +  &
1605            &        zCdn_form_ice * ( zfmi * at_i_b(ji,jj) + zfmw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) )
1606
1607         ! Heat Transfert Coefficients (Eq. 49)
1608         pch(ji,jj) = zChn_skin_ice *   zfhi +  &
1609            &        zChn_form_ice * ( zfhi * at_i_b(ji,jj) + zfhw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) )
1610         !
1611      END_2D
1612      CALL lbc_lnk_multi( 'sbcblk', pcd, 'T',  1.0_wp, pch, 'T', 1.0_wp )
1613      !
1614   END SUBROUTINE Cdn10_Lupkes2015
1615
1616#endif
1617
1618   !!======================================================================
1619END MODULE sbcblk
Note: See TracBrowser for help on using the repository browser.