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 @ 11405

Last change on this file since 11405 was 11405, checked in by andmirek, 5 years ago

ticket #2195: read weights for blk using XIOS

  • Property svn:keywords set to Id
File size: 77.4 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 (http://aerobulk.sourceforge.net/)
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   !!----------------------------------------------------------------------
21
22   !!----------------------------------------------------------------------
23   !!   sbc_blk_init  : initialisation of the chosen bulk formulation as ocean surface boundary condition
24   !!   sbc_blk       : bulk formulation as ocean surface boundary condition
25   !!   blk_oce       : computes momentum, heat and freshwater fluxes over ocean
26   !!   rho_air       : density of (moist) air (depends on T_air, q_air and SLP
27   !!   cp_air        : specific heat of (moist) air (depends spec. hum. q_air)
28   !!   q_sat         : saturation humidity as a function of SLP and temperature
29   !!   L_vap         : latent heat of vaporization of water as a function of temperature
30   !!             sea-ice case only :
31   !!   blk_ice_tau   : provide the air-ice stress
32   !!   blk_ice_flx   : provide the heat and mass fluxes at air-ice interface
33   !!   blk_ice_qcn   : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
34   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
35   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag
36   !!----------------------------------------------------------------------
37   USE oce            ! ocean dynamics and tracers
38   USE dom_oce        ! ocean space and time domain
39   USE phycst         ! physical constants
40   USE fldread        ! read input fields
41   USE sbc_oce        ! Surface boundary condition: ocean fields
42   USE cyclone        ! Cyclone 10m wind form trac of cyclone centres
43   USE sbcdcy         ! surface boundary condition: diurnal cycle
44   USE sbcwave , ONLY :   cdn_wave ! wave module
45   USE sbc_ice        ! Surface boundary condition: ice fields
46   USE lib_fortran    ! to use key_nosignedzero
47#if defined key_si3
48   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, t_su, rn_cnd_s, hfx_err_dif
49   USE icethd_dh      ! for CALL ice_thd_snwblow
50#endif
51   USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009)
52   USE sbcblk_algo_coare    ! => turb_coare    : COAREv3.0 (Fairall et al. 2003)
53   USE sbcblk_algo_coare3p5 ! => turb_coare3p5 : COAREv3.5 (Edson et al. 2013)
54   USE sbcblk_algo_ecmwf    ! => turb_ecmwf    : ECMWF (IFS cycle 31)
55   !
56   USE iom            ! I/O manager library
57   USE in_out_manager ! I/O manager
58   USE lib_mpp        ! distribued memory computing library
59   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
60   USE prtctl         ! Print control
61
62   IMPLICIT NONE
63   PRIVATE
64
65   PUBLIC   sbc_blk_init  ! called in sbcmod
66   PUBLIC   sbc_blk       ! called in sbcmod
67#if defined key_si3
68   PUBLIC   blk_ice_tau   ! routine called in icesbc
69   PUBLIC   blk_ice_flx   ! routine called in icesbc
70   PUBLIC   blk_ice_qcn   ! routine called in icesbc
71#endif 
72
73!!Lolo: should ultimately be moved in the module with all physical constants ?
74!!gm  : In principle, yes.
75   REAL(wp), PARAMETER ::   Cp_dry = 1005.0       !: Specic heat of dry air, constant pressure      [J/K/kg]
76   REAL(wp), PARAMETER ::   Cp_vap = 1860.0       !: Specic heat of water vapor, constant pressure  [J/K/kg]
77   REAL(wp), PARAMETER ::   R_dry = 287.05_wp     !: Specific gas constant for dry air              [J/K/kg]
78   REAL(wp), PARAMETER ::   R_vap = 461.495_wp    !: Specific gas constant for water vapor          [J/K/kg]
79   REAL(wp), PARAMETER ::   reps0 = R_dry/R_vap   !: ratio of gas constant for dry air and water vapor => ~ 0.622
80   REAL(wp), PARAMETER ::   rctv0 = R_vap/R_dry   !: for virtual temperature (== (1-eps)/eps) => ~ 0.608
81
82   INTEGER , PARAMETER ::   jpfld   =10           ! maximum number of files to read
83   INTEGER , PARAMETER ::   jp_wndi = 1           ! index of 10m wind velocity (i-component) (m/s)    at T-point
84   INTEGER , PARAMETER ::   jp_wndj = 2           ! index of 10m wind velocity (j-component) (m/s)    at T-point
85   INTEGER , PARAMETER ::   jp_tair = 3           ! index of 10m air temperature             (Kelvin)
86   INTEGER , PARAMETER ::   jp_humi = 4           ! index of specific humidity               ( % )
87   INTEGER , PARAMETER ::   jp_qsr  = 5           ! index of solar heat                      (W/m2)
88   INTEGER , PARAMETER ::   jp_qlw  = 6           ! index of Long wave                       (W/m2)
89   INTEGER , PARAMETER ::   jp_prec = 7           ! index of total precipitation (rain+snow) (Kg/m2/s)
90   INTEGER , PARAMETER ::   jp_snow = 8           ! index of snow (solid prcipitation)       (kg/m2/s)
91   INTEGER , PARAMETER ::   jp_slp  = 9           ! index of sea level pressure              (Pa)
92   INTEGER , PARAMETER ::   jp_tdif =10           ! index of tau diff associated to HF tau   (N/m2)   at T-point
93
94   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input fields (file informations, fields read)
95
96   !                                             !!! Bulk parameters
97   REAL(wp), PARAMETER ::   cpa    = 1000.5         ! specific heat of air (only used for ice fluxes now...)
98   REAL(wp), PARAMETER ::   Ls     =    2.839e6     ! latent heat of sublimation
99   REAL(wp), PARAMETER ::   Stef   =    5.67e-8     ! Stefan Boltzmann constant
100   REAL(wp), PARAMETER ::   Cd_ice =    1.4e-3      ! transfer coefficient over ice
101   REAL(wp), PARAMETER ::   albo   =    0.066       ! ocean albedo assumed to be constant
102   !
103   !                           !!* Namelist namsbc_blk : bulk parameters
104   LOGICAL  ::   ln_NCAR        ! "NCAR"      algorithm   (Large and Yeager 2008)
105   LOGICAL  ::   ln_COARE_3p0   ! "COARE 3.0" algorithm   (Fairall et al. 2003)
106   LOGICAL  ::   ln_COARE_3p5   ! "COARE 3.5" algorithm   (Edson et al. 2013)
107   LOGICAL  ::   ln_ECMWF       ! "ECMWF"     algorithm   (IFS cycle 31)
108   !
109   LOGICAL  ::   ln_taudif      ! logical flag to use the "mean of stress module - module of mean stress" data
110   REAL(wp) ::   rn_pfac        ! multiplication factor for precipitation
111   REAL(wp) ::   rn_efac        ! multiplication factor for evaporation
112   REAL(wp) ::   rn_vfac        ! multiplication factor for ice/ocean velocity in the calculation of wind stress
113   REAL(wp) ::   rn_zqt         ! z(q,t) : height of humidity and temperature measurements
114   REAL(wp) ::   rn_zu          ! z(u)   : height of wind measurements
115!!gm ref namelist initialize it so remove the setting to false below
116   LOGICAL  ::   ln_Cd_L12 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2012)
117   LOGICAL  ::   ln_Cd_L15 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2015)
118   !
119   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_atm                    ! transfer coefficient for momentum      (tau)
120   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ch_atm                    ! transfer coefficient for sensible heat (Q_sens)
121   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ce_atm                    ! tansfert coefficient for evaporation   (Q_lat)
122   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_zu                      ! air temperature at wind speed height (needed by Lupkes 2015 bulk scheme)
123   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   q_zu                      ! air spec. hum.  at wind speed height (needed by Lupkes 2015 bulk scheme)
124   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cdn_oce, chn_oce, cen_oce ! needed by Lupkes 2015 bulk scheme
125
126   INTEGER  ::   nblk           ! choice of the bulk algorithm
127   !                            ! associated indices:
128   INTEGER, PARAMETER ::   np_NCAR      = 1   ! "NCAR" algorithm        (Large and Yeager 2008)
129   INTEGER, PARAMETER ::   np_COARE_3p0 = 2   ! "COARE 3.0" algorithm   (Fairall et al. 2003)
130   INTEGER, PARAMETER ::   np_COARE_3p5 = 3   ! "COARE 3.5" algorithm   (Edson et al. 2013)
131   INTEGER, PARAMETER ::   np_ECMWF     = 4   ! "ECMWF" algorithm       (IFS cycle 31)
132
133   !! * Substitutions
134#  include "vectopt_loop_substitute.h90"
135   !!----------------------------------------------------------------------
136   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
137   !! $Id$
138   !! Software governed by the CeCILL license (see ./LICENSE)
139   !!----------------------------------------------------------------------
140CONTAINS
141
142   INTEGER FUNCTION sbc_blk_alloc()
143      !!-------------------------------------------------------------------
144      !!             ***  ROUTINE sbc_blk_alloc ***
145      !!-------------------------------------------------------------------
146      ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), &
147         &      cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc )
148      !
149      CALL mpp_sum ( 'sbcblk', sbc_blk_alloc )
150      IF( sbc_blk_alloc /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_alloc: failed to allocate arrays' )
151   END FUNCTION sbc_blk_alloc
152
153
154   SUBROUTINE sbc_blk_init
155      !!---------------------------------------------------------------------
156      !!                    ***  ROUTINE sbc_blk_init  ***
157      !!
158      !! ** Purpose :   choose and initialize a bulk formulae formulation
159      !!
160      !! ** Method  :
161      !!
162      !!----------------------------------------------------------------------
163      INTEGER  ::   ifpr, jfld            ! dummy loop indice and argument
164      INTEGER  ::   ios, ierror, ioptio   ! Local integer
165      !!
166      CHARACTER(len=100)            ::   cn_dir                ! Root directory for location of atmospheric forcing files
167      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                 ! array of namelist informations on the fields to read
168      TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_humi, sn_qsr       ! informations about the fields to be read
169      TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow      !       "                        "
170      TYPE(FLD_N) ::   sn_slp , sn_tdif                        !       "                        "
171      NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw ,                &   ! input fields
172         &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif,                &
173         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF,             &   ! bulk algorithm
174         &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         & 
175         &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15
176      !!---------------------------------------------------------------------
177      !
178      !                                      ! allocate sbc_blk_core array
179      IF( sbc_blk_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' )
180      !
181      !                             !** read bulk namelist 
182      REWIND( numnam_ref )                !* Namelist namsbc_blk in reference namelist : bulk parameters
183      READ  ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901)
184901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_blk in reference namelist', lwp )
185      !
186      REWIND( numnam_cfg )                !* Namelist namsbc_blk in configuration namelist : bulk parameters
187      READ  ( numnam_cfg, namsbc_blk, IOSTAT = ios, ERR = 902 )
188902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_blk in configuration namelist', lwp )
189      !
190      IF(lwm) WRITE( numond, namsbc_blk )
191      !
192      !                             !** initialization of the chosen bulk formulae (+ check)
193      !                                   !* select the bulk chosen in the namelist and check the choice
194                                                               ioptio = 0
195      IF( ln_NCAR      ) THEN   ;   nblk =  np_NCAR        ;   ioptio = ioptio + 1   ;   ENDIF
196      IF( ln_COARE_3p0 ) THEN   ;   nblk =  np_COARE_3p0   ;   ioptio = ioptio + 1   ;   ENDIF
197      IF( ln_COARE_3p5 ) THEN   ;   nblk =  np_COARE_3p5   ;   ioptio = ioptio + 1   ;   ENDIF
198      IF( ln_ECMWF     ) THEN   ;   nblk =  np_ECMWF       ;   ioptio = ioptio + 1   ;   ENDIF
199      !
200      IF( ioptio /= 1 )   CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' )
201      !
202      IF( ln_dm2dc ) THEN                 !* check: diurnal cycle on Qsr
203         IF( sn_qsr%nfreqh /= 24 )   CALL ctl_stop( 'sbc_blk_init: ln_dm2dc=T only with daily short-wave input' )
204         IF( sn_qsr%ln_tint ) THEN
205            CALL ctl_warn( 'sbc_blk_init: ln_dm2dc=T daily qsr time interpolation done by sbcdcy module',   &
206               &           '              ==> We force time interpolation = .false. for qsr' )
207            sn_qsr%ln_tint = .false.
208         ENDIF
209      ENDIF
210      !                                   !* set the bulk structure
211      !                                      !- store namelist information in an array
212      slf_i(jp_wndi) = sn_wndi   ;   slf_i(jp_wndj) = sn_wndj
213      slf_i(jp_qsr ) = sn_qsr    ;   slf_i(jp_qlw ) = sn_qlw
214      slf_i(jp_tair) = sn_tair   ;   slf_i(jp_humi) = sn_humi
215      slf_i(jp_prec) = sn_prec   ;   slf_i(jp_snow) = sn_snow
216      slf_i(jp_slp)  = sn_slp    ;   slf_i(jp_tdif) = sn_tdif
217      !
218      lhftau = ln_taudif                     !- add an extra field if HF stress is used
219      jfld = jpfld - COUNT( (/.NOT.lhftau/) )
220      !
221      !                                      !- allocate the bulk structure
222      ALLOCATE( sf(jfld), STAT=ierror )
223      IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_init: unable to allocate sf structure' )
224      DO ifpr= 1, jfld
225         ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
226         IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
227         IF( slf_i(ifpr)%nfreqh > 0. .AND. MOD( 3600. * slf_i(ifpr)%nfreqh , REAL(nn_fsbc) * rdt) /= 0. )   &
228            &  CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rdt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.',   &
229            &                 '               This is not ideal. You should consider changing either rdt or nn_fsbc value...' )
230
231      END DO
232      !                                      !- fill the bulk structure with namelist informations
233      CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' )
234      !
235      IF ( ln_wave ) THEN
236      !Activated wave module but neither drag nor stokes drift activated
237         IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) )   THEN
238            CALL ctl_stop( 'STOP',  'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' )
239      !drag coefficient read from wave model definable only with mfs bulk formulae and core
240         ELSEIF (ln_cdgw .AND. .NOT. ln_NCAR )       THEN       
241             CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae')
242         ELSEIF (ln_stcor .AND. .NOT. ln_sdw)                             THEN
243             CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')
244         ENDIF
245      ELSE
246      IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor )                & 
247         &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    &
248         &                  'with drag coefficient (ln_cdgw =T) '  ,                        &
249         &                  'or Stokes Drift (ln_sdw=T) ' ,                                 &
250         &                  'or ocean stress modification due to waves (ln_tauwoc=T) ',      & 
251         &                  'or Stokes-Coriolis term (ln_stcori=T)'  )
252      ENDIF 
253      !
254      !           
255      IF(lwp) THEN                     !** Control print
256         !
257         WRITE(numout,*)                  !* namelist
258         WRITE(numout,*) '   Namelist namsbc_blk (other than data information):'
259         WRITE(numout,*) '      "NCAR"      algorithm   (Large and Yeager 2008)     ln_NCAR      = ', ln_NCAR
260         WRITE(numout,*) '      "COARE 3.0" algorithm   (Fairall et al. 2003)       ln_COARE_3p0 = ', ln_COARE_3p0
261         WRITE(numout,*) '      "COARE 3.5" algorithm   (Edson et al. 2013)         ln_COARE_3p5 = ', ln_COARE_3p0
262         WRITE(numout,*) '      "ECMWF"     algorithm   (IFS cycle 31)              ln_ECMWF     = ', ln_ECMWF
263         WRITE(numout,*) '      add High freq.contribution to the stress module     ln_taudif    = ', ln_taudif
264         WRITE(numout,*) '      Air temperature and humidity reference height (m)   rn_zqt       = ', rn_zqt
265         WRITE(numout,*) '      Wind vector reference height (m)                    rn_zu        = ', rn_zu
266         WRITE(numout,*) '      factor applied on precipitation (total & snow)      rn_pfac      = ', rn_pfac
267         WRITE(numout,*) '      factor applied on evaporation                       rn_efac      = ', rn_efac
268         WRITE(numout,*) '      factor applied on ocean/ice velocity                rn_vfac      = ', rn_vfac
269         WRITE(numout,*) '         (form absolute (=0) to relative winds(=1))'
270         WRITE(numout,*) '      use ice-atm drag from Lupkes2012                    ln_Cd_L12    = ', ln_Cd_L12
271         WRITE(numout,*) '      use ice-atm drag from Lupkes2015                    ln_Cd_L15    = ', ln_Cd_L15
272         !
273         WRITE(numout,*)
274         SELECT CASE( nblk )              !* Print the choice of bulk algorithm
275         CASE( np_NCAR      )   ;   WRITE(numout,*) '   ==>>>   "NCAR" algorithm        (Large and Yeager 2008)'
276         CASE( np_COARE_3p0 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.0" algorithm   (Fairall et al. 2003)'
277         CASE( np_COARE_3p5 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.5" algorithm   (Edson et al. 2013)'
278         CASE( np_ECMWF     )   ;   WRITE(numout,*) '   ==>>>   "ECMWF" algorithm       (IFS cycle 31)'
279         END SELECT
280         !
281      ENDIF
282
283      IF(lxios_blkw) THEN
284         CALL iom_blk_wgt_init( cbwxios_context )
285         CALL iom_swap( cxios_context )
286      ENDIF
287     
288      !
289   END SUBROUTINE sbc_blk_init
290
291
292   SUBROUTINE sbc_blk( kt )
293      !!---------------------------------------------------------------------
294      !!                    ***  ROUTINE sbc_blk  ***
295      !!
296      !! ** Purpose :   provide at each time step the surface ocean fluxes
297      !!              (momentum, heat, freshwater and runoff)
298      !!
299      !! ** Method  : (1) READ each fluxes in NetCDF files:
300      !!      the 10m wind velocity (i-component) (m/s)    at T-point
301      !!      the 10m wind velocity (j-component) (m/s)    at T-point
302      !!      the 10m or 2m specific humidity     ( % )
303      !!      the solar heat                      (W/m2)
304      !!      the Long wave                       (W/m2)
305      !!      the 10m or 2m air temperature       (Kelvin)
306      !!      the total precipitation (rain+snow) (Kg/m2/s)
307      !!      the snow (solid prcipitation)       (kg/m2/s)
308      !!      the tau diff associated to HF tau   (N/m2)   at T-point   (ln_taudif=T)
309      !!              (2) CALL blk_oce
310      !!
311      !!      C A U T I O N : never mask the surface stress fields
312      !!                      the stress is assumed to be in the (i,j) mesh referential
313      !!
314      !! ** Action  :   defined at each time-step at the air-sea interface
315      !!              - utau, vtau  i- and j-component of the wind stress
316      !!              - taum        wind stress module at T-point
317      !!              - wndm        wind speed  module at T-point over free ocean or leads in presence of sea-ice
318      !!              - qns, qsr    non-solar and solar heat fluxes
319      !!              - emp         upward mass flux (evapo. - precip.)
320      !!              - sfx         salt flux due to freezing/melting (non-zero only if ice is present)
321      !!
322      !! ** References :   Large & Yeager, 2004 / Large & Yeager, 2008
323      !!                   Brodeau et al. Ocean Modelling 2010
324      !!----------------------------------------------------------------------
325      INTEGER, INTENT(in) ::   kt   ! ocean time step
326      !!---------------------------------------------------------------------
327      !
328      CALL fld_read( kt, nn_fsbc, sf )             ! input fields provided at the current time-step
329      !
330      !                                            ! compute the surface ocean fluxes using bulk formulea
331      IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce( kt, sf, sst_m, ssu_m, ssv_m )
332
333#if defined key_cice
334      IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN
335         qlw_ice(:,:,1)   = sf(jp_qlw )%fnow(:,:,1)
336         IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )
337         ELSE                ; qsr_ice(:,:,1) =          sf(jp_qsr)%fnow(:,:,1) 
338         ENDIF
339         tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)
340         qatm_ice(:,:)    = sf(jp_humi)%fnow(:,:,1)
341         tprecip(:,:)     = sf(jp_prec)%fnow(:,:,1) * rn_pfac
342         sprecip(:,:)     = sf(jp_snow)%fnow(:,:,1) * rn_pfac
343         wndi_ice(:,:)    = sf(jp_wndi)%fnow(:,:,1)
344         wndj_ice(:,:)    = sf(jp_wndj)%fnow(:,:,1)
345      ENDIF
346#endif
347      !
348   END SUBROUTINE sbc_blk
349
350
351   SUBROUTINE blk_oce( kt, sf, pst, pu, pv )
352      !!---------------------------------------------------------------------
353      !!                     ***  ROUTINE blk_oce  ***
354      !!
355      !! ** Purpose :   provide the momentum, heat and freshwater fluxes at
356      !!      the ocean surface at each time step
357      !!
358      !! ** Method  :   bulk formulea for the ocean using atmospheric
359      !!      fields read in sbc_read
360      !!
361      !! ** Outputs : - utau    : i-component of the stress at U-point  (N/m2)
362      !!              - vtau    : j-component of the stress at V-point  (N/m2)
363      !!              - taum    : Wind stress module at T-point         (N/m2)
364      !!              - wndm    : Wind speed module at T-point          (m/s)
365      !!              - qsr     : Solar heat flux over the ocean        (W/m2)
366      !!              - qns     : Non Solar heat flux over the ocean    (W/m2)
367      !!              - emp     : evaporation minus precipitation       (kg/m2/s)
368      !!
369      !!  ** Nota  :   sf has to be a dummy argument for AGRIF on NEC
370      !!---------------------------------------------------------------------
371      INTEGER  , INTENT(in   )                 ::   kt    ! time step index
372      TYPE(fld), INTENT(inout), DIMENSION(:)   ::   sf    ! input data
373      REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pst   ! surface temperature                      [Celcius]
374      REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pu    ! surface current at U-point (i-component) [m/s]
375      REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pv    ! surface current at V-point (j-component) [m/s]
376      !
377      INTEGER  ::   ji, jj               ! dummy loop indices
378      REAL(wp) ::   zztmp                ! local variable
379      REAL(wp), DIMENSION(jpi,jpj) ::   zwnd_i, zwnd_j    ! wind speed components at T-point
380      REAL(wp), DIMENSION(jpi,jpj) ::   zsq               ! specific humidity at pst
381      REAL(wp), DIMENSION(jpi,jpj) ::   zqlw, zqsb        ! long wave and sensible heat fluxes
382      REAL(wp), DIMENSION(jpi,jpj) ::   zqla, zevap       ! latent heat fluxes and evaporation
383      REAL(wp), DIMENSION(jpi,jpj) ::   zst               ! surface temperature in Kelvin
384      REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s]
385      REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K]
386      REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa             ! density of air   [kg/m^3]
387      !!---------------------------------------------------------------------
388      !
389      ! local scalars ( place there for vector optimisation purposes)
390      zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K)
391
392      ! ----------------------------------------------------------------------------- !
393      !      0   Wind components and module at T-point relative to the moving ocean   !
394      ! ----------------------------------------------------------------------------- !
395
396      ! ... components ( U10m - U_oce ) at T-point (unmasked)
397!!gm    move zwnd_i (_j) set to zero  inside the key_cyclone ???
398      zwnd_i(:,:) = 0._wp
399      zwnd_j(:,:) = 0._wp
400#if defined key_cyclone
401      CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012)
402      DO jj = 2, jpjm1
403         DO ji = fs_2, fs_jpim1   ! vect. opt.
404            sf(jp_wndi)%fnow(ji,jj,1) = sf(jp_wndi)%fnow(ji,jj,1) + zwnd_i(ji,jj)
405            sf(jp_wndj)%fnow(ji,jj,1) = sf(jp_wndj)%fnow(ji,jj,1) + zwnd_j(ji,jj)
406         END DO
407      END DO
408#endif
409      DO jj = 2, jpjm1
410         DO ji = fs_2, fs_jpim1   ! vect. opt.
411            zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  )
412            zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  )
413         END DO
414      END DO
415      CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. )
416      ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked)
417      wndm(:,:) = SQRT(  zwnd_i(:,:) * zwnd_i(:,:)   &
418         &             + zwnd_j(:,:) * zwnd_j(:,:)  ) * tmask(:,:,1)
419
420      ! ----------------------------------------------------------------------------- !
421      !      I   Radiative FLUXES                                                     !
422      ! ----------------------------------------------------------------------------- !
423
424      ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle                    ! Short Wave
425      zztmp = 1. - albo
426      IF( ln_dm2dc ) THEN   ;   qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)
427      ELSE                  ;   qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1)
428      ENDIF
429
430      zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave
431
432      ! ----------------------------------------------------------------------------- !
433      !     II    Turbulent FLUXES                                                    !
434      ! ----------------------------------------------------------------------------- !
435
436      ! ... specific humidity at SST and IST tmask(
437      zsq(:,:) = 0.98 * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) )
438      !!
439      !! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate
440      !!    (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2
441      !!    (since reanalysis products provide T at z, not theta !)
442      ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1) ) * rn_zqt
443
444      SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point
445      !
446      CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! NCAR-COREv2
447         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
448      CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.0
449         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
450      CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.5
451         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
452      CASE( np_ECMWF     )   ;   CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! ECMWF
453         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
454      CASE DEFAULT
455         CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' )
456      END SELECT
457
458      !                          ! Compute true air density :
459      IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN     ! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...)
460         zrhoa(:,:) = rho_air( t_zu(:,:)              , q_zu(:,:)              , sf(jp_slp)%fnow(:,:,1) )
461      ELSE                                      ! At zt:
462         zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
463      END IF
464
465!!      CALL iom_put( "Cd_oce", Cd_atm)  ! output value of pure ocean-atm. transfer coef.
466!!      CALL iom_put( "Ch_oce", Ch_atm)  ! output value of pure ocean-atm. transfer coef.
467
468      DO jj = 1, jpj             ! tau module, i and j component
469         DO ji = 1, jpi
470            zztmp = zrhoa(ji,jj)  * zU_zu(ji,jj) * Cd_atm(ji,jj)   ! using bulk wind speed
471            taum  (ji,jj) = zztmp * wndm  (ji,jj)
472            zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
473            zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
474         END DO
475      END DO
476
477      !                          ! add the HF tau contribution to the wind stress module
478      IF( lhftau )   taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1)
479
480      CALL iom_put( "taum_oce", taum )   ! output wind stress module
481
482      ! ... utau, vtau at U- and V_points, resp.
483      !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
484      !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
485      DO jj = 1, jpjm1
486         DO ji = 1, fs_jpim1
487            utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) ) &
488               &          * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))
489            vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) ) &
490               &          * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1))
491         END DO
492      END DO
493      CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. )
494
495      !  Turbulent fluxes over ocean
496      ! -----------------------------
497
498      ! zqla used as temporary array, for rho*U (common term of bulk formulae):
499      zqla(:,:) = zrhoa(:,:) * zU_zu(:,:) * tmask(:,:,1)
500
501      IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN
502         !! q_air and t_air are given at 10m (wind reference height)
503         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed
504         zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:)             )   ! Sensible Heat, using bulk wind speed
505      ELSE
506         !! q_air and t_air are not given at 10m (wind reference height)
507         ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!!
508         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed
509         zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) )   ! Sensible Heat, using bulk wind speed
510      ENDIF
511
512      zqla(:,:) = L_vap(zst(:,:)) * zevap(:,:)     ! Latent Heat flux
513
514
515      IF(ln_ctl) THEN
516         CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce: zqla   : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce  : ' )
517         CALL prt_ctl( tab2d_1=zqsb  , clinfo1=' blk_oce: zqsb   : ', tab2d_2=Ch_atm , clinfo2=' Ch_oce  : ' )
518         CALL prt_ctl( tab2d_1=zqlw  , clinfo1=' blk_oce: zqlw   : ', tab2d_2=qsr, clinfo2=' qsr : ' )
519         CALL prt_ctl( tab2d_1=zsq   , clinfo1=' blk_oce: zsq    : ', tab2d_2=zst, clinfo2=' zst : ' )
520         CALL prt_ctl( tab2d_1=utau  , clinfo1=' blk_oce: utau   : ', mask1=umask,   &
521            &          tab2d_2=vtau  , clinfo2=           ' vtau : ', mask2=vmask )
522         CALL prt_ctl( tab2d_1=wndm  , clinfo1=' blk_oce: wndm   : ')
523         CALL prt_ctl( tab2d_1=zst   , clinfo1=' blk_oce: zst    : ')
524      ENDIF
525
526      ! ----------------------------------------------------------------------------- !
527      !     III    Total FLUXES                                                       !
528      ! ----------------------------------------------------------------------------- !
529      !
530      emp (:,:) = (  zevap(:,:)                                          &   ! mass flux (evap. - precip.)
531         &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1)
532      !
533      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar
534         &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * rLfus                        &   ! remove latent melting heat for solid precip
535         &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST
536         &     + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac  &   ! add liquid precip heat content at Tair
537         &     * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp                          &
538         &     + sf(jp_snow)%fnow(:,:,1) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow)
539         &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi
540      qns(:,:) = qns(:,:) * tmask(:,:,1)
541      !
542#if defined key_si3
543      qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                ! non solar without emp (only needed by SI3)
544      qsr_oce(:,:) = qsr(:,:)
545#endif
546      !
547      IF ( nn_ice == 0 ) THEN
548         CALL iom_put( "qlw_oce" ,   zqlw )                 ! output downward longwave heat over the ocean
549         CALL iom_put( "qsb_oce" , - zqsb )                 ! output downward sensible heat over the ocean
550         CALL iom_put( "qla_oce" , - zqla )                 ! output downward latent   heat over the ocean
551         CALL iom_put( "qemp_oce",   qns-zqlw+zqsb+zqla )   ! output downward heat content of E-P over the ocean
552         CALL iom_put( "qns_oce" ,   qns  )                 ! output downward non solar heat over the ocean
553         CALL iom_put( "qsr_oce" ,   qsr  )                 ! output downward solar heat over the ocean
554         CALL iom_put( "qt_oce"  ,   qns+qsr )              ! output total downward heat over the ocean
555         tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s]
556         sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s]
557         CALL iom_put( 'snowpre', sprecip )                 ! Snow
558         CALL iom_put( 'precip' , tprecip )                 ! Total precipitation
559      ENDIF
560      !
561      IF(ln_ctl) THEN
562         CALL prt_ctl(tab2d_1=zqsb , clinfo1=' blk_oce: zqsb   : ', tab2d_2=zqlw , clinfo2=' zqlw  : ')
563         CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce: zqla   : ', tab2d_2=qsr  , clinfo2=' qsr   : ')
564         CALL prt_ctl(tab2d_1=pst  , clinfo1=' blk_oce: pst    : ', tab2d_2=emp  , clinfo2=' emp   : ')
565         CALL prt_ctl(tab2d_1=utau , clinfo1=' blk_oce: utau   : ', mask1=umask,   &
566            &         tab2d_2=vtau , clinfo2=              ' vtau  : ' , mask2=vmask )
567      ENDIF
568      !
569   END SUBROUTINE blk_oce
570
571
572
573   FUNCTION rho_air( ptak, pqa, pslp )
574      !!-------------------------------------------------------------------------------
575      !!                           ***  FUNCTION rho_air  ***
576      !!
577      !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere
578      !!
579      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
580      !!-------------------------------------------------------------------------------
581      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K]
582      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg]
583      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp      ! pressure in                [Pa]
584      REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air   ! density of moist air   [kg/m^3]
585      !!-------------------------------------------------------------------------------
586      !
587      rho_air = pslp / (  R_dry*ptak * ( 1._wp + rctv0*pqa )  )
588      !
589   END FUNCTION rho_air
590
591
592   FUNCTION cp_air( pqa )
593      !!-------------------------------------------------------------------------------
594      !!                           ***  FUNCTION cp_air  ***
595      !!
596      !! ** Purpose : Compute specific heat (Cp) of moist air
597      !!
598      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
599      !!-------------------------------------------------------------------------------
600      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg]
601      REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air   ! specific heat of moist air   [J/K/kg]
602      !!-------------------------------------------------------------------------------
603      !
604      Cp_air = Cp_dry + Cp_vap * pqa
605      !
606   END FUNCTION cp_air
607
608
609   FUNCTION q_sat( ptak, pslp )
610      !!----------------------------------------------------------------------------------
611      !!                           ***  FUNCTION q_sat  ***
612      !!
613      !! ** Purpose : Specific humidity at saturation in [kg/kg]
614      !!              Based on accurate estimate of "e_sat"
615      !!              aka saturation water vapor (Goff, 1957)
616      !!
617      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
618      !!----------------------------------------------------------------------------------
619      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak    ! air temperature                       [K]
620      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp    ! sea level atmospheric pressure       [Pa]
621      REAL(wp), DIMENSION(jpi,jpj)             ::   q_sat   ! Specific humidity at saturation   [kg/kg]
622      !
623      INTEGER  ::   ji, jj         ! dummy loop indices
624      REAL(wp) ::   ze_sat, ztmp   ! local scalar
625      !!----------------------------------------------------------------------------------
626      !
627      DO jj = 1, jpj
628         DO ji = 1, jpi
629            !
630            ztmp = rt0 / ptak(ji,jj)
631            !
632            ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957)
633            ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0)        &
634               &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) )  &
635               &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614  )
636               !
637            q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wp*pslp(ji,jj) - (1._wp - reps0)*ze_sat )   ! 0.01 because SLP is in [Pa]
638            !
639         END DO
640      END DO
641      !
642   END FUNCTION q_sat
643
644
645   FUNCTION gamma_moist( ptak, pqa )
646      !!----------------------------------------------------------------------------------
647      !!                           ***  FUNCTION gamma_moist  ***
648      !!
649      !! ** Purpose : Compute the moist adiabatic lapse-rate.
650      !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
651      !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html
652      !!
653      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
654      !!----------------------------------------------------------------------------------
655      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak          ! air temperature       [K]
656      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! specific humidity [kg/kg]
657      REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist   ! moist adiabatic lapse-rate
658      !
659      INTEGER  ::   ji, jj         ! dummy loop indices
660      REAL(wp) :: zrv, ziRT        ! local scalar
661      !!----------------------------------------------------------------------------------
662      !
663      DO jj = 1, jpj
664         DO ji = 1, jpi
665            zrv = pqa(ji,jj) / (1. - pqa(ji,jj))
666            ziRT = 1. / (R_dry*ptak(ji,jj))    ! 1/RT
667            gamma_moist(ji,jj) = grav * ( 1. + rLevap*zrv*ziRT ) / ( Cp_dry + rLevap*rLevap*zrv*reps0*ziRT/ptak(ji,jj) )
668         END DO
669      END DO
670      !
671   END FUNCTION gamma_moist
672
673
674   FUNCTION L_vap( psst )
675      !!---------------------------------------------------------------------------------
676      !!                           ***  FUNCTION L_vap  ***
677      !!
678      !! ** Purpose : Compute the latent heat of vaporization of water from temperature
679      !!
680      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
681      !!----------------------------------------------------------------------------------
682      REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! latent heat of vaporization   [J/kg]
683      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K]
684      !!----------------------------------------------------------------------------------
685      !
686      L_vap = (  2.501 - 0.00237 * ( psst(:,:) - rt0)  ) * 1.e6
687      !
688   END FUNCTION L_vap
689
690   SUBROUTINE iom_blk_wgt_init( cdname, ld_tmppatch ) 
691#if defined key_iomput
692      use xios
693#endif
694      !!----------------------------------------------------------------------
695      !!                     ***  ROUTINE   ***
696      !!
697      !! ** Purpose :   initialize context for reading weights for surface
698      !!                forcing
699      !!
700      !!----------------------------------------------------------------------
701      CHARACTER(len=*),           INTENT(in)  :: cdname
702      LOGICAL         , OPTIONAL, INTENT(in)  :: ld_tmppatch
703#if defined key_iomput
704      !
705      INTEGER, PARAMETER :: lcname = 26
706      CHARACTER(len=lc) :: clname
707      INTEGER           :: ji, jkmin
708      LOGICAL ::   ll_tmppatch = .TRUE.    !: seb: patch before we remove periodicity
709      INTEGER ::   nldi_save, nlei_save    !:      and close boundaries in output files
710      INTEGER ::   nldj_save, nlej_save    !:
711      LOGICAL ::   ll_global = .FALSE.     !: do we have variable on model grid
712      CHARACTER(len=lc), DIMENSION( jpts) :: cg_name(jpts)
713      CHARACTER(len=1), DIMENSION(lcname) :: cname
714      CHARACTER(len=lc)               :: cfname ! file name without .nc
715      TYPE(xios_duration)               :: dtime    = xios_duration(0, 0, 0, 0, 0, 0), &
716                                           outp_frq = xios_duration(0, 0, 0, 0, 0, 0)
717      TYPE(xios_domaingroup)            :: domaingroup_hdl
718      TYPE(xios_domain)                 :: domain_hdl
719      TYPE(xios_axisgroup)              :: axisgroup_hdl
720      TYPE(xios_axis)                   :: axis_hdl
721      TYPE(xios_scalar)                 :: scalar_hdl
722      TYPE(xios_scalargroup)            :: scalargroup_hdl
723      TYPE(xios_file)                   :: file_hdl
724      TYPE(xios_filegroup)              :: filegroup_hdl
725      TYPE(xios_field)                  :: field_hdl
726      INTEGER                           :: jf, ni, nj, ipos, jfld
727      INTEGER, DIMENSION(3)             :: ndims         ! size of each dimension. this is surface
728                                                         ! forcing, not more than 3: (x, y, t)
729      INTEGER                           :: ndim          ! number if dimensions
730      INTEGER                           :: inum          ! file id
731      INTEGER                           :: ivid          ! varable id in netcdf file
732      CHARACTER (len=5)                 :: aname         !name of the field in wght file
733      INTEGER                           :: numwgt        !
734      LOGICAL                           :: lfound
735      INTEGER                           :: jg, id, jn
736     
737
738      cname(:)= (/'a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', &
739                  'k', 'l', 'm', 'n', 'o', 'p', 'q', 'r', 's', 't', &
740                  'u', 'v', 'w', 'x', 'y', 'z'/)
741      !!----------------------------------------------------------------------
742      !
743      ! seb: patch before we remove periodicity and close boundaries in output files
744      IF ( ll_tmppatch ) THEN
745         nldi_save = nldi   ;   nlei_save = nlei
746         nldj_save = nldj   ;   nlej_save = nlej
747         IF( nimpp           ==      1 ) nldi = 1
748         IF( nimpp + jpi - 1 == jpiglo ) nlei = jpi
749         IF( njmpp           ==      1 ) nldj = 1
750         IF( njmpp + jpj - 1 == jpjglo ) nlej = jpj
751      ENDIF
752      !
753      jfld = SIZE(sf)
754
755      IF(jfld > lcname) THEN
756      STOP 'lcname number is too small'
757      ENDIF
758
759      clname = cdname
760      IF( TRIM(Agrif_CFixed()) /= '0' )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(cdname)
761      CALL xios_context_initialize(TRIM(clname), mpi_comm_oce)
762      CALL iom_swap( cdname )
763      ! Calendar type is now defined in xml file
764      SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
765      CASE ( 1)   ; CALL xios_define_calendar( TYPE = "Gregorian", time_origin = xios_date(1948, 01, 01, 00, 00, 00), &
766          &                                    start_date = xios_date(nyear, nmonth, nday,00,00,00) )
767      CASE ( 0)   ; CALL xios_define_calendar( TYPE = "NoLeap"   , time_origin = xios_date(1948, 01, 01, 00, 00, 00), &
768          &                                    start_date = xios_date(nyear, nmonth, nday,00,00,00) )
769      CASE (30)   ; CALL xios_define_calendar( TYPE = "D360"     , time_origin = xios_date(1948, 01, 01, 00, 00, 00), &
770          &                                    start_date = xios_date(nyear, nmonth, nday, 00, 00, 00) )
771      END SELECT
772
773
774      dtime%month = 1
775      CALL xios_set_timestep( dtime )
776
777      CALL xios_get_handle("domain_definition",domaingroup_hdl)
778      CALL xios_add_child(domaingroup_hdl, domain_hdl, "grid_global")
779      ni = nlei-nldi+1
780      nj = nlej-nldj+1
781      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)
782      CALL xios_set_domain_attr("grid_global", data_dim=2, data_ibegin = 1-nldi, data_ni = jpi, data_jbegin = 1-nldj, data_nj = jpj)
783      CALL xios_set_domain_attr("grid_global", type='curvilinear')
784
785      CALL xios_get_handle("file_definition", filegroup_hdl )
786
787      DO jf = 1, jfld
788!do we have weights file defined already?
789         lfound = .FALSE.
790         DO jg = 1, jf-1
791            IF(TRIM(sf(jf)%wgtname) == TRIM(sf(jg)%wgtname)) THEN
792               IF(lwp) write(numout, *) TRIM(sf(jf)%wgtname),' already defined, skipping '
793               lfound = .TRUE.
794               exit
795            ENDIF
796         ENDDO
797         IF(lfound) cycle
798!add next weights file
799         ipos = index(sf(jf)%wgtname,'.nc')
800         cfname(1:lc) = " "
801         IF(ipos > 0) THEN
802          cfname(1:ipos-1) = sf(jf)%wgtname(1:ipos-1)
803         ELSE
804          cfname(1:lc) = sf(jf)%wgtname(1:lc)
805         ENDIF
806
807         CALL xios_add_child(filegroup_hdl, file_hdl, cname(jf))
808         CALL xios_set_file_attr( cname(jf), name=TRIM(cfname),           &
809                type="one_file", par_access="collective", enabled=.TRUE., &
810                output_freq=xios_timestep, mode="read")
811
812         CALL iom_open ( sf(jf)%wgtname, inum )
813         id = iom_varid(inum, 'src05', ldstop=.FALSE.)
814         IF(id <= 0) THEN
815             numwgt = 4
816         ELSE
817             numwgt = 16
818         ENDIF
819         CALL iom_close ( inum )
820!define variables in weights file. Because all files have the same names src01,
821!dst01, ... add filename to the name to uniquely identify variables
822         aname = ' '
823         DO jn = 1, numwgt
824            WRITE(aname,'(a3,i2.2)') 'src',jn
825            CALL xios_add_child(file_hdl, field_hdl, TRIM(aname)//TRIM(sf(jf)%wgtname))
826            CALL xios_set_attr (field_hdl, enabled = .TRUE., &
827                          name = TRIM(aname), domain_ref="grid_global", &
828                          operation = "instant")
829            CALL xios_add_child(file_hdl, field_hdl, TRIM(aname)//TRIM(sf(jf)%wgtname))
830            CALL xios_set_attr (field_hdl, enabled = .TRUE., &
831                          name = TRIM(aname), domain_ref="grid_global", &
832                          operation = "instant") 
833         ENDDO 
834      ENDDO   
835
836
837      CALL xios_close_context_definition()
838      CALL xios_update_calendar( 0 )
839
840      IF ( ll_tmppatch ) THEN
841         nldi = nldi_save   ;   nlei = nlei_save
842         nldj = nldj_save   ;   nlej = nlej_save
843      ENDIF
844#endif
845      !
846   END SUBROUTINE iom_blk_wgt_init
847
848
849#if defined key_si3
850   !!----------------------------------------------------------------------
851   !!   'key_si3'                                       SI3 sea-ice model
852   !!----------------------------------------------------------------------
853   !!   blk_ice_tau : provide the air-ice stress
854   !!   blk_ice_flx : provide the heat and mass fluxes at air-ice interface
855   !!   blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
856   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
857   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag
858   !!----------------------------------------------------------------------
859
860   SUBROUTINE blk_ice_tau
861      !!---------------------------------------------------------------------
862      !!                     ***  ROUTINE blk_ice_tau  ***
863      !!
864      !! ** Purpose :   provide the surface boundary condition over sea-ice
865      !!
866      !! ** Method  :   compute momentum using bulk formulation
867      !!                formulea, ice variables and read atmospheric fields.
868      !!                NB: ice drag coefficient is assumed to be a constant
869      !!---------------------------------------------------------------------
870      INTEGER  ::   ji, jj    ! dummy loop indices
871      REAL(wp) ::   zwndi_f , zwndj_f, zwnorm_f   ! relative wind module and components at F-point
872      REAL(wp) ::   zwndi_t , zwndj_t             ! relative wind components at T-point
873      REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa     ! transfer coefficient for momentum      (tau)
874      !!---------------------------------------------------------------------
875      !
876      ! set transfer coefficients to default sea-ice values
877      Cd_atm(:,:) = Cd_ice
878      Ch_atm(:,:) = Cd_ice
879      Ce_atm(:,:) = Cd_ice
880
881      wndm_ice(:,:) = 0._wp      !!gm brutal....
882
883      ! ------------------------------------------------------------ !
884      !    Wind module relative to the moving ice ( U10m - U_ice )   !
885      ! ------------------------------------------------------------ !
886      ! C-grid ice dynamics :   U & V-points (same as ocean)
887      DO jj = 2, jpjm1
888         DO ji = fs_2, fs_jpim1   ! vect. opt.
889            zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  )
890            zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  )
891            wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
892         END DO
893      END DO
894      CALL lbc_lnk( 'sbcblk', wndm_ice, 'T',  1. )
895      !
896      ! Make ice-atm. drag dependent on ice concentration
897      IF    ( ln_Cd_L12 ) THEN   ! calculate new drag from Lupkes(2012) equations
898         CALL Cdn10_Lupkes2012( Cd_atm )
899         Ch_atm(:,:) = Cd_atm(:,:)       ! momentum and heat transfer coef. are considered identical
900      ELSEIF( ln_Cd_L15 ) THEN   ! calculate new drag from Lupkes(2015) equations
901         CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm ) 
902      ENDIF
903
904!!      CALL iom_put( "Cd_ice", Cd_atm)  ! output value of pure ice-atm. transfer coef.
905!!      CALL iom_put( "Ch_ice", Ch_atm)  ! output value of pure ice-atm. transfer coef.
906
907      ! local scalars ( place there for vector optimisation purposes)
908      ! Computing density of air! Way denser that 1.2 over sea-ice !!!
909      zrhoa (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1))
910
911      !!gm brutal....
912      utau_ice  (:,:) = 0._wp
913      vtau_ice  (:,:) = 0._wp
914      !!gm end
915
916      ! ------------------------------------------------------------ !
917      !    Wind stress relative to the moving ice ( U10m - U_ice )   !
918      ! ------------------------------------------------------------ !
919      ! C-grid ice dynamics :   U & V-points (same as ocean)
920      DO jj = 2, jpjm1
921         DO ji = fs_2, fs_jpim1   ! vect. opt.
922            utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )            &
923               &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) )
924            vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )            &
925               &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) )
926         END DO
927      END DO
928      CALL lbc_lnk_multi( 'sbcblk', utau_ice, 'U', -1., vtau_ice, 'V', -1. )
929      !
930      !
931      IF(ln_ctl) THEN
932         CALL prt_ctl(tab2d_1=utau_ice  , clinfo1=' blk_ice: utau_ice : ', tab2d_2=vtau_ice  , clinfo2=' vtau_ice : ')
933         CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice: wndm_ice : ')
934      ENDIF
935      !
936   END SUBROUTINE blk_ice_tau
937
938
939   SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb )
940      !!---------------------------------------------------------------------
941      !!                     ***  ROUTINE blk_ice_flx  ***
942      !!
943      !! ** Purpose :   provide the heat and mass fluxes at air-ice interface
944      !!
945      !! ** Method  :   compute heat and freshwater exchanged
946      !!                between atmosphere and sea-ice using bulk formulation
947      !!                formulea, ice variables and read atmmospheric fields.
948      !!
949      !! caution : the net upward water flux has with mm/day unit
950      !!---------------------------------------------------------------------
951      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu   ! sea ice surface temperature
952      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness
953      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness
954      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb   ! ice albedo (all skies)
955      !!
956      INTEGER  ::   ji, jj, jl               ! dummy loop indices
957      REAL(wp) ::   zst3                     ! local variable
958      REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      -
959      REAL(wp) ::   zztmp, z1_rLsub           !   -      -
960      REAL(wp) ::   zfr1, zfr2               ! local variables
961      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_st         ! inverse of surface temperature
962      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qlw         ! long wave heat flux over ice
963      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qsb         ! sensible  heat flux over ice
964      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_dqlw        ! long wave heat sensitivity over ice
965      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_dqsb        ! sensible  heat sensitivity over ice
966      REAL(wp), DIMENSION(jpi,jpj)     ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (SI3)
967      REAL(wp), DIMENSION(jpi,jpj)     ::   zrhoa
968      !!---------------------------------------------------------------------
969      !
970      zcoef_dqlw = 4.0 * 0.95 * Stef             ! local scalars
971      zcoef_dqla = -Ls * 11637800. * (-5897.8)
972      !
973      zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
974      !
975      zztmp = 1. / ( 1. - albo )
976      WHERE( ptsu(:,:,:) /= 0._wp )   ;   z1_st(:,:,:) = 1._wp / ptsu(:,:,:)
977      ELSEWHERE                       ;   z1_st(:,:,:) = 0._wp
978      END WHERE
979      !                                     ! ========================== !
980      DO jl = 1, jpl                        !  Loop over ice categories  !
981         !                                  ! ========================== !
982         DO jj = 1 , jpj
983            DO ji = 1, jpi
984               ! ----------------------------!
985               !      I   Radiative FLUXES   !
986               ! ----------------------------!
987               zst3 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl) * ptsu(ji,jj,jl)
988               ! Short Wave (sw)
989               qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
990               ! Long  Wave (lw)
991               z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
992               ! lw sensitivity
993               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
994
995               ! ----------------------------!
996               !     II    Turbulent FLUXES  !
997               ! ----------------------------!
998
999               ! ... turbulent heat fluxes with Ch_atm recalculated in blk_ice_tau
1000               ! Sensible Heat
1001               z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1))
1002               ! Latent Heat
1003               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls  * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  &
1004                  &                ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) )
1005               ! Latent heat sensitivity for ice (Dqla/Dt)
1006               IF( qla_ice(ji,jj,jl) > 0._wp ) THEN
1007                  dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  &
1008                     &                 z1_st(ji,jj,jl)*z1_st(ji,jj,jl) * EXP(-5897.8 * z1_st(ji,jj,jl))
1009               ELSE
1010                  dqla_ice(ji,jj,jl) = 0._wp
1011               ENDIF
1012
1013               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice)
1014               z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj)
1015
1016               ! ----------------------------!
1017               !     III    Total FLUXES     !
1018               ! ----------------------------!
1019               ! Downward Non Solar flux
1020               qns_ice (ji,jj,jl) =     z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl)
1021               ! Total non solar heat flux sensitivity for ice
1022               dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) )
1023            END DO
1024            !
1025         END DO
1026         !
1027      END DO
1028      !
1029      tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1)  ! total precipitation [kg/m2/s]
1030      sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1)  ! solid precipitation [kg/m2/s]
1031      CALL iom_put( 'snowpre', sprecip )                    ! Snow precipitation
1032      CALL iom_put( 'precip' , tprecip )                    ! Total precipitation
1033
1034      ! --- evaporation --- !
1035      z1_rLsub = 1._wp / rLsub
1036      evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_rLsub    ! sublimation
1037      devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_rLsub    ! d(sublimation)/dT
1038      zevap    (:,:)   = rn_efac * ( emp(:,:) + tprecip(:,:) )   ! evaporation over ocean
1039
1040      ! --- evaporation minus precipitation --- !
1041      zsnw(:,:) = 0._wp
1042      CALL ice_thd_snwblow( (1.-at_i_b(:,:)), zsnw )  ! snow distribution over ice after wind blowing
1043      emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
1044      emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
1045      emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:)
1046
1047      ! --- heat flux associated with emp --- !
1048      qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp                  & ! evap at sst
1049         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip at Tair
1050         &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip at min(Tair,Tsnow)
1051         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
1052      qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only)
1053         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
1054
1055      ! --- total solar and non solar fluxes --- !
1056      qns_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 )  &
1057         &           + qemp_ice(:,:) + qemp_oce(:,:)
1058      qsr_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
1059
1060      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
1061      qprec_ice(:,:) = rhos * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
1062
1063      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) ---
1064      DO jl = 1, jpl
1065         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * rcpi * tmask(:,:,1) )
1066         !                         ! But we do not have Tice => consider it at 0degC => evap=0
1067      END DO
1068
1069      ! --- shortwave radiation transmitted below the surface (W/m2, see Grenfell Maykut 77) --- !
1070      zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            ! transmission when hi>10cm
1071      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            ! zfr2 such that zfr1 + zfr2 to equal 1
1072      !
1073      WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm 
1074         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) )
1075      ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (zfr1) when hi>10cm
1076         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1
1077      ELSEWHERE                                                         ! zero when hs>0
1078         qtr_ice_top(:,:,:) = 0._wp 
1079      END WHERE
1080      !
1081      IF(ln_ctl) THEN
1082         CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb    : ', kdim=jpl)
1083         CALL prt_ctl(tab3d_1=z_qlw   , clinfo1=' blk_ice: z_qlw    : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl)
1084         CALL prt_ctl(tab3d_1=z_dqsb  , clinfo1=' blk_ice: z_dqsb   : ', tab3d_2=z_dqlw  , clinfo2=' z_dqlw   : ', kdim=jpl)
1085         CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice: dqns_ice : ', tab3d_2=qsr_ice , clinfo2=' qsr_ice  : ', kdim=jpl)
1086         CALL prt_ctl(tab3d_1=ptsu    , clinfo1=' blk_ice: ptsu     : ', tab3d_2=qns_ice , clinfo2=' qns_ice  : ', kdim=jpl)
1087         CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip  : ', tab2d_2=sprecip , clinfo2=' sprecip  : ')
1088      ENDIF
1089      !
1090   END SUBROUTINE blk_ice_flx
1091   
1092
1093   SUBROUTINE blk_ice_qcn( ld_virtual_itd, ptsu, ptb, phs, phi )
1094      !!---------------------------------------------------------------------
1095      !!                     ***  ROUTINE blk_ice_qcn  ***
1096      !!
1097      !! ** Purpose :   Compute surface temperature and snow/ice conduction flux
1098      !!                to force sea ice / snow thermodynamics
1099      !!                in the case conduction flux is emulated
1100      !!               
1101      !! ** Method  :   compute surface energy balance assuming neglecting heat storage
1102      !!                following the 0-layer Semtner (1976) approach
1103      !!
1104      !! ** Outputs : - ptsu    : sea-ice / snow surface temperature (K)
1105      !!              - qcn_ice : surface inner conduction flux (W/m2)
1106      !!
1107      !!---------------------------------------------------------------------
1108      LOGICAL                   , INTENT(in   ) ::   ld_virtual_itd  ! single-category option
1109      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptsu            ! sea ice / snow surface temperature
1110      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   ptb             ! sea ice base temperature
1111      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   phs             ! snow thickness
1112      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   phi             ! sea ice thickness
1113      !
1114      INTEGER , PARAMETER ::   nit = 10                  ! number of iterations
1115      REAL(wp), PARAMETER ::   zepsilon = 0.1_wp         ! characteristic thickness for enhanced conduction
1116      !
1117      INTEGER  ::   ji, jj, jl           ! dummy loop indices
1118      INTEGER  ::   iter                 ! local integer
1119      REAL(wp) ::   zfac, zfac2, zfac3   ! local scalars
1120      REAL(wp) ::   zkeff_h, ztsu, ztsu0 !
1121      REAL(wp) ::   zqc, zqnet           !
1122      REAL(wp) ::   zhe, zqa0            !
1123      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zgfac   ! enhanced conduction factor
1124      !!---------------------------------------------------------------------
1125     
1126      ! -------------------------------------!
1127      !      I   Enhanced conduction factor  !
1128      ! -------------------------------------!
1129      ! Emulates the enhancement of conduction by unresolved thin ice (ld_virtual_itd = T)
1130      ! Fichefet and Morales Maqueda, JGR 1997
1131      !
1132      zgfac(:,:,:) = 1._wp
1133     
1134      IF( ld_virtual_itd ) THEN
1135         !
1136         zfac  = 1._wp /  ( rn_cnd_s + rcnd_i )
1137         zfac2 = EXP(1._wp) * 0.5_wp * zepsilon
1138         zfac3 = 2._wp / zepsilon
1139         !   
1140         DO jl = 1, jpl               
1141            DO jj = 1 , jpj
1142               DO ji = 1, jpi
1143                  zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac                            ! Effective thickness
1144                  IF( zhe >=  zfac2 )   zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor
1145               END DO
1146            END DO
1147         END DO
1148         !     
1149      ENDIF
1150     
1151      ! -------------------------------------------------------------!
1152      !      II   Surface temperature and conduction flux            !
1153      ! -------------------------------------------------------------!
1154      !
1155      zfac = rcnd_i * rn_cnd_s
1156      !
1157      DO jl = 1, jpl
1158         DO jj = 1 , jpj
1159            DO ji = 1, jpi
1160               !                   
1161               zkeff_h = zfac * zgfac(ji,jj,jl) / &                                    ! Effective conductivity of the snow-ice system divided by thickness
1162                  &      ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) )
1163               ztsu    = ptsu(ji,jj,jl)                                                ! Store current iteration temperature
1164               ztsu0   = ptsu(ji,jj,jl)                                                ! Store initial surface temperature
1165               zqa0    = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux
1166               !
1167               DO iter = 1, nit     ! --- Iterative loop
1168                  zqc   = zkeff_h * ( ztsu - ptb(ji,jj) )                              ! Conduction heat flux through snow-ice system (>0 downwards)
1169                  zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc  ! Surface energy budget
1170                  ztsu  = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h )              ! Temperature update
1171               END DO
1172               !
1173               ptsu   (ji,jj,jl) = MIN( rt0, ztsu )
1174               qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) )
1175               qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 )
1176               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) )  &
1177                             &   * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) )
1178
1179               ! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- !
1180               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) 
1181
1182            END DO
1183         END DO
1184         !
1185      END DO 
1186      !     
1187   END SUBROUTINE blk_ice_qcn
1188   
1189
1190   SUBROUTINE Cdn10_Lupkes2012( Cd )
1191      !!----------------------------------------------------------------------
1192      !!                      ***  ROUTINE  Cdn10_Lupkes2012  ***
1193      !!
1194      !! ** Purpose :    Recompute the neutral air-ice drag referenced at 10m
1195      !!                 to make it dependent on edges at leads, melt ponds and flows.
1196      !!                 After some approximations, this can be resumed to a dependency
1197      !!                 on ice concentration.
1198      !!               
1199      !! ** Method :     The parameterization is taken from Lupkes et al. (2012) eq.(50)
1200      !!                 with the highest level of approximation: level4, eq.(59)
1201      !!                 The generic drag over a cell partly covered by ice can be re-written as follows:
1202      !!
1203      !!                 Cd = Cdw * (1-A) + Cdi * A + Ce * (1-A)**(nu+1/(10*beta)) * A**mu
1204      !!
1205      !!                    Ce = 2.23e-3       , as suggested by Lupkes (eq. 59)
1206      !!                    nu = mu = beta = 1 , as suggested by Lupkes (eq. 59)
1207      !!                    A is the concentration of ice minus melt ponds (if any)
1208      !!
1209      !!                 This new drag has a parabolic shape (as a function of A) starting at
1210      !!                 Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5
1211      !!                 and going down to Cdi(say 1.4e-3) for A=1
1212      !!
1213      !!                 It is theoretically applicable to all ice conditions (not only MIZ)
1214      !!                 => see Lupkes et al (2013)
1215      !!
1216      !! ** References : Lupkes et al. JGR 2012 (theory)
1217      !!                 Lupkes et al. GRL 2013 (application to GCM)
1218      !!
1219      !!----------------------------------------------------------------------
1220      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd
1221      REAL(wp), PARAMETER ::   zCe   = 2.23e-03_wp
1222      REAL(wp), PARAMETER ::   znu   = 1._wp
1223      REAL(wp), PARAMETER ::   zmu   = 1._wp
1224      REAL(wp), PARAMETER ::   zbeta = 1._wp
1225      REAL(wp)            ::   zcoef
1226      !!----------------------------------------------------------------------
1227      zcoef = znu + 1._wp / ( 10._wp * zbeta )
1228
1229      ! generic drag over a cell partly covered by ice
1230      !!Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) +  &                        ! pure ocean drag
1231      !!   &      Cd_ice      *           at_i_b(:,:)   +  &                        ! pure ice drag
1232      !!   &      zCe         * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu   ! change due to sea-ice morphology
1233
1234      ! ice-atm drag
1235      Cd(:,:) = Cd_ice +  &                                                         ! pure ice drag
1236         &      zCe    * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp)  ! change due to sea-ice morphology
1237     
1238   END SUBROUTINE Cdn10_Lupkes2012
1239
1240
1241   SUBROUTINE Cdn10_Lupkes2015( Cd, Ch )
1242      !!----------------------------------------------------------------------
1243      !!                      ***  ROUTINE  Cdn10_Lupkes2015  ***
1244      !!
1245      !! ** pUrpose :    Alternative turbulent transfert coefficients formulation
1246      !!                 between sea-ice and atmosphere with distinct momentum
1247      !!                 and heat coefficients depending on sea-ice concentration
1248      !!                 and atmospheric stability (no meltponds effect for now).
1249      !!               
1250      !! ** Method :     The parameterization is adapted from Lupkes et al. (2015)
1251      !!                 and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme,
1252      !!                 it considers specific skin and form drags (Andreas et al. 2010)
1253      !!                 to compute neutral transfert coefficients for both heat and
1254      !!                 momemtum fluxes. Atmospheric stability effect on transfert
1255      !!                 coefficient is also taken into account following Louis (1979).
1256      !!
1257      !! ** References : Lupkes et al. JGR 2015 (theory)
1258      !!                 Lupkes et al. ECHAM6 documentation 2015 (implementation)
1259      !!
1260      !!----------------------------------------------------------------------
1261      !
1262      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd
1263      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Ch
1264      REAL(wp), DIMENSION(jpi,jpj)            ::   ztm_su, zst, zqo_sat, zqi_sat
1265      !
1266      ! ECHAM6 constants
1267      REAL(wp), PARAMETER ::   z0_skin_ice  = 0.69e-3_wp  ! Eq. 43 [m]
1268      REAL(wp), PARAMETER ::   z0_form_ice  = 0.57e-3_wp  ! Eq. 42 [m]
1269      REAL(wp), PARAMETER ::   z0_ice       = 1.00e-3_wp  ! Eq. 15 [m]
1270      REAL(wp), PARAMETER ::   zce10        = 2.80e-3_wp  ! Eq. 41
1271      REAL(wp), PARAMETER ::   zbeta        = 1.1_wp      ! Eq. 41
1272      REAL(wp), PARAMETER ::   zc           = 5._wp       ! Eq. 13
1273      REAL(wp), PARAMETER ::   zc2          = zc * zc
1274      REAL(wp), PARAMETER ::   zam          = 2. * zc     ! Eq. 14
1275      REAL(wp), PARAMETER ::   zah          = 3. * zc     ! Eq. 30
1276      REAL(wp), PARAMETER ::   z1_alpha     = 1._wp / 0.2_wp  ! Eq. 51
1277      REAL(wp), PARAMETER ::   z1_alphaf    = z1_alpha    ! Eq. 56
1278      REAL(wp), PARAMETER ::   zbetah       = 1.e-3_wp    ! Eq. 26
1279      REAL(wp), PARAMETER ::   zgamma       = 1.25_wp     ! Eq. 26
1280      REAL(wp), PARAMETER ::   z1_gamma     = 1._wp / zgamma
1281      REAL(wp), PARAMETER ::   r1_3         = 1._wp / 3._wp
1282      !
1283      INTEGER  ::   ji, jj         ! dummy loop indices
1284      REAL(wp) ::   zthetav_os, zthetav_is, zthetav_zu
1285      REAL(wp) ::   zrib_o, zrib_i
1286      REAL(wp) ::   zCdn_skin_ice, zCdn_form_ice, zCdn_ice
1287      REAL(wp) ::   zChn_skin_ice, zChn_form_ice
1288      REAL(wp) ::   z0w, z0i, zfmi, zfmw, zfhi, zfhw
1289      REAL(wp) ::   zCdn_form_tmp
1290      !!----------------------------------------------------------------------
1291
1292      ! mean temperature
1293      WHERE( at_i_b(:,:) > 1.e-20 )   ;   ztm_su(:,:) = SUM( t_su(:,:,:) * a_i_b(:,:,:) , dim=3 ) / at_i_b(:,:)
1294      ELSEWHERE                       ;   ztm_su(:,:) = rt0
1295      ENDWHERE
1296     
1297      ! Momentum Neutral Transfert Coefficients (should be a constant)
1298      zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2   ! Eq. 40
1299      zCdn_skin_ice = ( vkarmn                                      / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2   ! Eq. 7
1300      zCdn_ice      = zCdn_skin_ice   ! Eq. 7 (cf Lupkes email for details)
1301      !zCdn_ice     = 1.89e-3         ! old ECHAM5 value (cf Eq. 32)
1302
1303      ! Heat Neutral Transfert Coefficients
1304      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 (cf Lupkes email for details)
1305     
1306      ! Atmospheric and Surface Variables
1307      zst(:,:)     = sst_m(:,:) + rt0                                        ! convert SST from Celcius to Kelvin
1308      zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:)   , sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ocean [kg/kg]
1309      zqi_sat(:,:) = 0.98_wp * q_sat( ztm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg]
1310      !
1311      DO jj = 2, jpjm1           ! reduced loop is necessary for reproducibility
1312         DO ji = fs_2, fs_jpim1
1313            ! Virtual potential temperature [K]
1314            zthetav_os = zst(ji,jj)    * ( 1._wp + rctv0 * zqo_sat(ji,jj) )   ! over ocean
1315            zthetav_is = ztm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) )   ! ocean ice
1316            zthetav_zu = t_zu (ji,jj)  * ( 1._wp + rctv0 * q_zu(ji,jj)    )   ! at zu
1317           
1318            ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead)
1319            zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj)     )**2   ! over ocean
1320            zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2   ! over ice
1321           
1322            ! Momentum and Heat Neutral Transfert Coefficients
1323            zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta  ! Eq. 40
1324            zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) )               ! Eq. 53
1325                       
1326            ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead)
1327            z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water
1328            z0i = z0_skin_ice                                             ! over ice (cf Lupkes email for details)
1329            IF( zrib_o <= 0._wp ) THEN
1330               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
1331               zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )**r1_3 / ( Chn_oce(ji,jj) * MAX(0.01, wndm(ji,jj)) )   &     ! Eq. 26
1332                  &             )**zgamma )**z1_gamma
1333            ELSE
1334               zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 12
1335               zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 28
1336            ENDIF
1337           
1338            IF( zrib_i <= 0._wp ) THEN
1339               zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq.  9
1340               zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq. 25
1341            ELSE
1342               zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 11
1343               zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 27
1344            ENDIF
1345           
1346            ! Momentum Transfert Coefficients (Eq. 38)
1347            Cd(ji,jj) = zCdn_skin_ice *   zfmi +  &
1348               &        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) )
1349           
1350            ! Heat Transfert Coefficients (Eq. 49)
1351            Ch(ji,jj) = zChn_skin_ice *   zfhi +  &
1352               &        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) )
1353            !
1354         END DO
1355      END DO
1356      CALL lbc_lnk_multi( 'sbcblk', Cd, 'T',  1., Ch, 'T', 1. )
1357      !
1358   END SUBROUTINE Cdn10_Lupkes2015
1359
1360#endif
1361
1362   !!======================================================================
1363END MODULE sbcblk
Note: See TracBrowser for help on using the repository browser.