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.
sbccpl.F90 in branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 8755

Last change on this file since 8755 was 8755, checked in by jcastill, 6 years ago

Further changes for ticket #1980
Receive the ocean wind stress components from a wave model, both in forced and coupled mode

File size: 151.9 KB
Line 
1MODULE sbccpl
2   !!======================================================================
3   !!                       ***  MODULE  sbccpl  ***
4   !! Surface Boundary Condition :  momentum, heat and freshwater fluxes in coupled mode
5   !!======================================================================
6   !! History :  2.0  ! 2007-06  (R. Redler, N. Keenlyside, W. Park) Original code split into flxmod & taumod
7   !!            3.0  ! 2008-02  (G. Madec, C Talandier)  surface module
8   !!            3.1  ! 2009_02  (G. Madec, S. Masson, E. Maisonave, A. Caubel) generic coupled interface
9   !!            3.4  ! 2011_11  (C. Harris) more flexibility + multi-category fields
10   !!----------------------------------------------------------------------
11   !!----------------------------------------------------------------------
12   !!   namsbc_cpl      : coupled formulation namlist
13   !!   sbc_cpl_init    : initialisation of the coupled exchanges
14   !!   sbc_cpl_rcv     : receive fields from the atmosphere over the ocean (ocean only)
15   !!                     receive stress from the atmosphere over the ocean (ocean-ice case)
16   !!   sbc_cpl_ice_tau : receive stress from the atmosphere over ice
17   !!   sbc_cpl_ice_flx : receive fluxes from the atmosphere over ice
18   !!   sbc_cpl_snd     : send     fields to the atmosphere
19   !!----------------------------------------------------------------------
20   USE dom_oce         ! ocean space and time domain
21   USE sbc_oce         ! Surface boundary condition: ocean fields
22   USE trc_oce         ! share SMS/Ocean variables
23   USE sbc_ice         ! Surface boundary condition: ice fields
24   USE sbcapr          ! Stochastic param. : ???
25   USE sbcdcy          ! surface boundary condition: diurnal cycle
26   USE sbcwave         ! surface boundary condition: waves
27   USE phycst          ! physical constants
28#if defined key_lim3
29   USE ice            ! ice variables
30#endif
31#if defined key_lim2
32   USE par_ice_2      ! ice parameters
33   USE ice_2          ! ice variables
34#endif
35   USE cpl_oasis3     ! OASIS3 coupling
36   USE geo2ocean      !
37   USE oce   , ONLY : tsn, un, vn, sshn, ub, vb, sshb, fraqsr_1lev
38   USE albedo         !
39   USE eosbn2         !
40   USE sbcrnf, ONLY : l_rnfcpl
41   USE sbcisf   , ONLY : l_isfcpl
42#if defined key_cice
43   USE ice_domain_size, only: ncat
44#endif
45#if defined key_lim3
46   USE limthd_dh      ! for CALL lim_thd_snwblow
47#endif
48   !
49   USE in_out_manager ! I/O manager
50   USE iom            ! NetCDF library
51   USE lib_mpp        ! distribued memory computing library
52   USE wrk_nemo       ! work arrays
53   USE timing         ! Timing
54   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
55
56   IMPLICIT NONE
57   PRIVATE
58
59   PUBLIC   sbc_cpl_init      ! routine called by sbcmod.F90
60   PUBLIC   sbc_cpl_rcv       ! routine called by sbc_ice_lim(_2).F90
61   PUBLIC   sbc_cpl_snd       ! routine called by step.F90
62   PUBLIC   sbc_cpl_ice_tau   ! routine called by sbc_ice_lim(_2).F90
63   PUBLIC   sbc_cpl_ice_flx   ! routine called by sbc_ice_lim(_2).F90
64   PUBLIC   sbc_cpl_alloc     ! routine called in sbcice_cice.F90
65
66   INTEGER, PARAMETER ::   jpr_otx1   =  1   ! 3 atmosphere-ocean stress components on grid 1
67   INTEGER, PARAMETER ::   jpr_oty1   =  2   !
68   INTEGER, PARAMETER ::   jpr_otz1   =  3   !
69   INTEGER, PARAMETER ::   jpr_otx2   =  4   ! 3 atmosphere-ocean stress components on grid 2
70   INTEGER, PARAMETER ::   jpr_oty2   =  5   !
71   INTEGER, PARAMETER ::   jpr_otz2   =  6   !
72   INTEGER, PARAMETER ::   jpr_itx1   =  7   ! 3 atmosphere-ice   stress components on grid 1
73   INTEGER, PARAMETER ::   jpr_ity1   =  8   !
74   INTEGER, PARAMETER ::   jpr_itz1   =  9   !
75   INTEGER, PARAMETER ::   jpr_itx2   = 10   ! 3 atmosphere-ice   stress components on grid 2
76   INTEGER, PARAMETER ::   jpr_ity2   = 11   !
77   INTEGER, PARAMETER ::   jpr_itz2   = 12   !
78   INTEGER, PARAMETER ::   jpr_qsroce = 13   ! Qsr above the ocean
79   INTEGER, PARAMETER ::   jpr_qsrice = 14   ! Qsr above the ice
80   INTEGER, PARAMETER ::   jpr_qsrmix = 15 
81   INTEGER, PARAMETER ::   jpr_qnsoce = 16   ! Qns above the ocean
82   INTEGER, PARAMETER ::   jpr_qnsice = 17   ! Qns above the ice
83   INTEGER, PARAMETER ::   jpr_qnsmix = 18
84   INTEGER, PARAMETER ::   jpr_rain   = 19   ! total liquid precipitation (rain)
85   INTEGER, PARAMETER ::   jpr_snow   = 20   ! solid precipitation over the ocean (snow)
86   INTEGER, PARAMETER ::   jpr_tevp   = 21   ! total evaporation
87   INTEGER, PARAMETER ::   jpr_ievp   = 22   ! solid evaporation (sublimation)
88   INTEGER, PARAMETER ::   jpr_sbpr   = 23   ! sublimation - liquid precipitation - solid precipitation
89   INTEGER, PARAMETER ::   jpr_semp   = 24   ! solid freshwater budget (sublimation - snow)
90   INTEGER, PARAMETER ::   jpr_oemp   = 25   ! ocean freshwater budget (evap - precip)
91   INTEGER, PARAMETER ::   jpr_w10m   = 26   ! 10m wind
92   INTEGER, PARAMETER ::   jpr_dqnsdt = 27   ! d(Q non solar)/d(temperature)
93   INTEGER, PARAMETER ::   jpr_rnf    = 28   ! runoffs
94   INTEGER, PARAMETER ::   jpr_cal    = 29   ! calving
95   INTEGER, PARAMETER ::   jpr_taum   = 30   ! wind stress module
96   INTEGER, PARAMETER ::   jpr_co2    = 31
97   INTEGER, PARAMETER ::   jpr_topm   = 32   ! topmeltn
98   INTEGER, PARAMETER ::   jpr_botm   = 33   ! botmeltn
99   INTEGER, PARAMETER ::   jpr_sflx   = 34   ! salt flux
100   INTEGER, PARAMETER ::   jpr_toce   = 35   ! ocean temperature
101   INTEGER, PARAMETER ::   jpr_soce   = 36   ! ocean salinity
102   INTEGER, PARAMETER ::   jpr_ocx1   = 37   ! ocean current on grid 1
103   INTEGER, PARAMETER ::   jpr_ocy1   = 38   !
104   INTEGER, PARAMETER ::   jpr_ssh    = 39   ! sea surface height
105   INTEGER, PARAMETER ::   jpr_fice   = 40   ! ice fraction         
106   INTEGER, PARAMETER ::   jpr_e3t1st = 41   ! first T level thickness
107   INTEGER, PARAMETER ::   jpr_fraqsr = 42   ! fraction of solar net radiation absorbed in the first ocean level
108   INTEGER, PARAMETER ::   jpr_mslp   = 43   ! mean sea level pressure
109   INTEGER, PARAMETER ::   jpr_hsig   = 44   ! Hsig
110   INTEGER, PARAMETER ::   jpr_phioc  = 45   ! Wave=>ocean energy flux
111   INTEGER, PARAMETER ::   jpr_sdrftx = 46   ! Stokes drift on grid 1
112   INTEGER, PARAMETER ::   jpr_sdrfty = 47   ! Stokes drift on grid 2
113   INTEGER, PARAMETER ::   jpr_wper   = 48   ! Mean wave period
114   INTEGER, PARAMETER ::   jpr_wnum   = 49   ! Mean wavenumber
115   INTEGER, PARAMETER ::   jpr_tauoc  = 50   ! Stress fraction adsorbed by waves
116   INTEGER, PARAMETER ::   jpr_wdrag  = 51   ! Neutral surface drag coefficient
117   INTEGER, PARAMETER ::   jpr_isf    = 52
118   INTEGER, PARAMETER ::   jpr_icb    = 53
119   INTEGER, PARAMETER ::   jpr_wfreq  = 54   ! Wave peak frequency
120   INTEGER, PARAMETER ::   jpr_tauwx  = 55   ! x component of the ocean stress from waves
121   INTEGER, PARAMETER ::   jpr_tauwy  = 56   ! y component of the ocean stress from waves
122
123   INTEGER, PARAMETER ::   jprcv      = 56   ! total number of fields received 
124
125   INTEGER, PARAMETER ::   jps_fice   =  1   ! ice fraction sent to the atmosphere
126   INTEGER, PARAMETER ::   jps_toce   =  2   ! ocean temperature
127   INTEGER, PARAMETER ::   jps_tice   =  3   ! ice   temperature
128   INTEGER, PARAMETER ::   jps_tmix   =  4   ! mixed temperature (ocean+ice)
129   INTEGER, PARAMETER ::   jps_albice =  5   ! ice   albedo
130   INTEGER, PARAMETER ::   jps_albmix =  6   ! mixed albedo
131   INTEGER, PARAMETER ::   jps_hice   =  7   ! ice  thickness
132   INTEGER, PARAMETER ::   jps_hsnw   =  8   ! snow thickness
133   INTEGER, PARAMETER ::   jps_ocx1   =  9   ! ocean current on grid 1
134   INTEGER, PARAMETER ::   jps_ocy1   = 10   !
135   INTEGER, PARAMETER ::   jps_ocz1   = 11   !
136   INTEGER, PARAMETER ::   jps_ivx1   = 12   ! ice   current on grid 1
137   INTEGER, PARAMETER ::   jps_ivy1   = 13   !
138   INTEGER, PARAMETER ::   jps_ivz1   = 14   !
139   INTEGER, PARAMETER ::   jps_co2    = 15
140   INTEGER, PARAMETER ::   jps_soce   = 16   ! ocean salinity
141   INTEGER, PARAMETER ::   jps_ssh    = 17   ! sea surface height
142   INTEGER, PARAMETER ::   jps_qsroce = 18   ! Qsr above the ocean
143   INTEGER, PARAMETER ::   jps_qnsoce = 19   ! Qns above the ocean
144   INTEGER, PARAMETER ::   jps_oemp   = 20   ! ocean freshwater budget (evap - precip)
145   INTEGER, PARAMETER ::   jps_sflx   = 21   ! salt flux
146   INTEGER, PARAMETER ::   jps_otx1   = 22   ! 2 atmosphere-ocean stress components on grid 1
147   INTEGER, PARAMETER ::   jps_oty1   = 23   !
148   INTEGER, PARAMETER ::   jps_rnf    = 24   ! runoffs
149   INTEGER, PARAMETER ::   jps_taum   = 25   ! wind stress module
150   INTEGER, PARAMETER ::   jps_fice2  = 26   ! ice fraction sent to OPA (by SAS when doing SAS-OPA coupling)
151   INTEGER, PARAMETER ::   jps_e3t1st = 27   ! first level depth (vvl)
152   INTEGER, PARAMETER ::   jps_fraqsr = 28   ! fraction of solar net radiation absorbed in the first ocean level
153   INTEGER, PARAMETER ::   jps_ficet  = 29   ! total ice fraction 
154   INTEGER, PARAMETER ::   jps_ocxw   = 30   ! currents on grid 1 
155   INTEGER, PARAMETER ::   jps_ocyw   = 31   ! currents on grid 2
156   INTEGER, PARAMETER ::   jps_wlev   = 32   ! water level
157   INTEGER, PARAMETER ::   jpsnd      = 32   ! total number of fields sent
158
159   !                                  !!** namelist namsbc_cpl **
160   TYPE ::   FLD_C                     !   
161      CHARACTER(len = 32) ::   cldes      ! desciption of the coupling strategy
162      CHARACTER(len = 32) ::   clcat      ! multiple ice categories strategy
163      CHARACTER(len = 32) ::   clvref     ! reference of vector ('spherical' or 'cartesian')
164      CHARACTER(len = 32) ::   clvor      ! orientation of vector fields ('eastward-northward' or 'local grid')
165      CHARACTER(len = 32) ::   clvgrd     ! grids on which is located the vector fields
166   END TYPE FLD_C
167   !                                   ! Send to the atmosphere 
168   TYPE(FLD_C) ::   sn_snd_temp, sn_snd_alb, sn_snd_thick, sn_snd_crt, sn_snd_co2                       
169   !                                   ! Received from the atmosphere
170   TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_tauw, sn_rcv_dqnsdt, sn_rcv_qsr, &
171                                                            sn_rcv_qns, sn_rcv_emp, sn_rcv_rnf
172   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_mslp, sn_rcv_icb, sn_rcv_isf                             
173   ! Send to waves
174   TYPE(FLD_C) ::   sn_snd_ifrac, sn_snd_crtw, sn_snd_wlev 
175   ! Received from waves
176   TYPE(FLD_C) ::   sn_rcv_hsig,sn_rcv_phioc,sn_rcv_sdrfx,sn_rcv_sdrfy,sn_rcv_wper,sn_rcv_wnum,sn_rcv_tauoc,sn_rcv_wdrag, &
177                    sn_rcv_wfreq
178   !                                   ! Other namelist parameters
179   INTEGER     ::   nn_cplmodel           ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
180   LOGICAL     ::   ln_usecplmask         !  use a coupling mask file to merge data received from several models
181                                         !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
182   TYPE ::   DYNARR     
183      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z3   
184   END TYPE DYNARR
185
186   TYPE( DYNARR ), SAVE, DIMENSION(jprcv) ::   frcv                     ! all fields recieved from the atmosphere
187
188   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   albedo_oce_mix    ! ocean albedo sent to atmosphere (mix clear/overcast sky)
189
190   REAL(wp) ::   rpref = 101000._wp   ! reference atmospheric pressure[N/m2]
191   REAL(wp) ::   r1_grau              ! = 1.e0 / (grav * rau0)
192
193   INTEGER , ALLOCATABLE, SAVE, DIMENSION(    :) ::   nrcvinfo           ! OASIS info argument
194
195   !! Substitution
196#  include "vectopt_loop_substitute.h90"
197   !!----------------------------------------------------------------------
198   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
199   !! $Id$
200   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
201   !!----------------------------------------------------------------------
202CONTAINS
203 
204   INTEGER FUNCTION sbc_cpl_alloc()
205      !!----------------------------------------------------------------------
206      !!             ***  FUNCTION sbc_cpl_alloc  ***
207      !!----------------------------------------------------------------------
208      INTEGER :: ierr(4)
209      !!----------------------------------------------------------------------
210      ierr(:) = 0
211      !
212      ALLOCATE( albedo_oce_mix(jpi,jpj), nrcvinfo(jprcv),  STAT=ierr(1) )
213     
214#if ! defined key_lim3 && ! defined key_lim2 && ! defined key_cice
215      ALLOCATE( a_i(jpi,jpj,1) , STAT=ierr(2) )  ! used in sbcice_if.F90 (done here as there is no sbc_ice_if_init)
216#endif
217      ALLOCATE( xcplmask(jpi,jpj,0:nn_cplmodel) , STAT=ierr(3) )
218      !
219      IF( .NOT. ln_apr_dyn ) ALLOCATE( ssh_ib(jpi,jpj), ssh_ibb(jpi,jpj), apr(jpi, jpj), STAT=ierr(4) ) 
220
221      sbc_cpl_alloc = MAXVAL( ierr )
222      IF( lk_mpp            )   CALL mpp_sum ( sbc_cpl_alloc )
223      IF( sbc_cpl_alloc > 0 )   CALL ctl_warn('sbc_cpl_alloc: allocation of arrays failed')
224      !
225   END FUNCTION sbc_cpl_alloc
226
227
228   SUBROUTINE sbc_cpl_init( k_ice )     
229      !!----------------------------------------------------------------------
230      !!             ***  ROUTINE sbc_cpl_init  ***
231      !!
232      !! ** Purpose :   Initialisation of send and received information from
233      !!                the atmospheric component
234      !!
235      !! ** Method  : * Read namsbc_cpl namelist
236      !!              * define the receive interface
237      !!              * define the send    interface
238      !!              * initialise the OASIS coupler
239      !!----------------------------------------------------------------------
240      INTEGER, INTENT(in) ::   k_ice   ! ice management in the sbc (=0/1/2/3)
241      !
242      INTEGER ::   jn          ! dummy loop index
243      INTEGER ::   ios, inum   ! Local integer
244      REAL(wp), POINTER, DIMENSION(:,:) ::   zacs, zaos
245      !!
246      NAMELIST/namsbc_cpl/  sn_snd_temp , sn_snd_alb  , sn_snd_thick , sn_snd_crt   , sn_snd_co2,      & 
247         &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau   , sn_rcv_dqnsdt, sn_rcv_qsr,      & 
248         &                  sn_snd_ifrac, sn_snd_crtw , sn_snd_wlev  , sn_rcv_hsig  , sn_rcv_phioc ,   & 
249         &                  sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper  , sn_rcv_wnum  , sn_rcv_tauoc ,   &
250         &                  sn_rcv_wdrag, sn_rcv_qns  , sn_rcv_emp   , sn_rcv_rnf   , sn_rcv_cal   ,   &
251         &                  sn_rcv_iceflx,sn_rcv_co2  , nn_cplmodel  , ln_usecplmask, sn_rcv_mslp  ,   &
252         &                  sn_rcv_icb , sn_rcv_isf   , sn_rcv_wfreq , sn_rcv_tauw
253
254      !!---------------------------------------------------------------------
255      !
256      IF( nn_timing == 1 )   CALL timing_start('sbc_cpl_init')
257      !
258      CALL wrk_alloc( jpi,jpj,   zacs, zaos )
259
260      ! ================================ !
261      !      Namelist informations       !
262      ! ================================ !
263      !
264      REWIND( numnam_ref )              ! Namelist namsbc_cpl in reference namelist : Variables for OASIS coupling
265      READ  ( numnam_ref, namsbc_cpl, IOSTAT = ios, ERR = 901)
266901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_cpl in reference namelist', lwp )
267      !
268      REWIND( numnam_cfg )              ! Namelist namsbc_cpl in configuration namelist : Variables for OASIS coupling
269      READ  ( numnam_cfg, namsbc_cpl, IOSTAT = ios, ERR = 902 )
270902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_cpl in configuration namelist', lwp )
271      IF(lwm) WRITE ( numond, namsbc_cpl )
272      !
273      IF(lwp) THEN                        ! control print
274         WRITE(numout,*)
275         WRITE(numout,*)'sbc_cpl_init : namsbc_cpl namelist '
276         WRITE(numout,*)'~~~~~~~~~~~~'
277      ENDIF
278      IF( lwp .AND. ln_cpl ) THEN                        ! control print
279         WRITE(numout,*)'  received fields (mutiple ice categogies)'
280         WRITE(numout,*)'      10m wind module                 = ', TRIM(sn_rcv_w10m%cldes  ), ' (', TRIM(sn_rcv_w10m%clcat  ), ')'
281         WRITE(numout,*)'      stress module                   = ', TRIM(sn_rcv_taumod%cldes), ' (', TRIM(sn_rcv_taumod%clcat), ')'
282         WRITE(numout,*)'      surface stress                  = ', TRIM(sn_rcv_tau%cldes   ), ' (', TRIM(sn_rcv_tau%clcat   ), ')'
283         WRITE(numout,*)'                     - referential    = ', sn_rcv_tau%clvref
284         WRITE(numout,*)'                     - orientation    = ', sn_rcv_tau%clvor
285         WRITE(numout,*)'                     - mesh           = ', sn_rcv_tau%clvgrd
286         WRITE(numout,*)'      non-solar heat flux sensitivity = ', TRIM(sn_rcv_dqnsdt%cldes), ' (', TRIM(sn_rcv_dqnsdt%clcat), ')'
287         WRITE(numout,*)'      solar heat flux                 = ', TRIM(sn_rcv_qsr%cldes   ), ' (', TRIM(sn_rcv_qsr%clcat   ), ')'
288         WRITE(numout,*)'      non-solar heat flux             = ', TRIM(sn_rcv_qns%cldes   ), ' (', TRIM(sn_rcv_qns%clcat   ), ')'
289         WRITE(numout,*)'      freshwater budget               = ', TRIM(sn_rcv_emp%cldes   ), ' (', TRIM(sn_rcv_emp%clcat   ), ')'
290         WRITE(numout,*)'      runoffs                         = ', TRIM(sn_rcv_rnf%cldes   ), ' (', TRIM(sn_rcv_rnf%clcat   ), ')'
291         WRITE(numout,*)'      calving                         = ', TRIM(sn_rcv_cal%cldes   ), ' (', TRIM(sn_rcv_cal%clcat   ), ')'
292         WRITE(numout,*)'      iceberg                         = ', TRIM(sn_rcv_icb%cldes   ), ' (', TRIM(sn_rcv_icb%clcat   ), ')'
293         WRITE(numout,*)'      ice shelf                       = ', TRIM(sn_rcv_isf%cldes   ), ' (', TRIM(sn_rcv_isf%clcat   ), ')'
294         WRITE(numout,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')'
295         WRITE(numout,*)'      atm co2                         = ', TRIM(sn_rcv_co2%cldes   ), ' (', TRIM(sn_rcv_co2%clcat   ), ')'
296         WRITE(numout,*)'      significant wave heigth         = ', TRIM(sn_rcv_hsig%cldes  ), ' (', TRIM(sn_rcv_hsig%clcat  ), ')' 
297         WRITE(numout,*)'      wave to oce energy flux         = ', TRIM(sn_rcv_phioc%cldes ), ' (', TRIM(sn_rcv_phioc%clcat ), ')' 
298         WRITE(numout,*)'      Surface Stokes drift grid u     = ', TRIM(sn_rcv_sdrfx%cldes ), ' (', TRIM(sn_rcv_sdrfx%clcat ), ')' 
299         WRITE(numout,*)'      Surface Stokes drift grid v     = ', TRIM(sn_rcv_sdrfy%cldes ), ' (', TRIM(sn_rcv_sdrfy%clcat ), ')' 
300         WRITE(numout,*)'      Mean wave period                = ', TRIM(sn_rcv_wper%cldes  ), ' (', TRIM(sn_rcv_wper%clcat  ), ')' 
301         WRITE(numout,*)'      Mean wave number                = ', TRIM(sn_rcv_wnum%cldes  ), ' (', TRIM(sn_rcv_wnum%clcat  ), ')' 
302         WRITE(numout,*)'      Wave peak frequency             = ', TRIM(sn_rcv_wfreq%cldes ), ' (', TRIM(sn_rcv_wfreq%clcat ), ')'
303         WRITE(numout,*)'      Stress frac adsorbed by waves   = ', TRIM(sn_rcv_tauoc%cldes ), ' (', TRIM(sn_rcv_tauoc%clcat ), ')' 
304         WRITE(numout,*)'      Stress components by waves      = ', TRIM(sn_rcv_tauw%cldes  ), ' (', TRIM(sn_rcv_tauw%clcat  ), ')'
305         WRITE(numout,*)'      Neutral surf drag coefficient   = ', TRIM(sn_rcv_wdrag%cldes ), ' (', TRIM(sn_rcv_wdrag%clcat ), ')' 
306         WRITE(numout,*)'  sent fields (multiple ice categories)'
307         WRITE(numout,*)'      surface temperature             = ', TRIM(sn_snd_temp%cldes  ), ' (', TRIM(sn_snd_temp%clcat  ), ')'
308         WRITE(numout,*)'      albedo                          = ', TRIM(sn_snd_alb%cldes   ), ' (', TRIM(sn_snd_alb%clcat   ), ')'
309         WRITE(numout,*)'      ice/snow thickness              = ', TRIM(sn_snd_thick%cldes ), ' (', TRIM(sn_snd_thick%clcat ), ')'
310         WRITE(numout,*)'      total ice fraction              = ', TRIM(sn_snd_ifrac%cldes ), ' (', TRIM(sn_snd_ifrac%clcat ), ')' 
311         WRITE(numout,*)'      surface current                 = ', TRIM(sn_snd_crt%cldes   ), ' (', TRIM(sn_snd_crt%clcat   ), ')'
312         WRITE(numout,*)'                      - referential   = ', sn_snd_crt%clvref 
313         WRITE(numout,*)'                      - orientation   = ', sn_snd_crt%clvor
314         WRITE(numout,*)'                      - mesh          = ', sn_snd_crt%clvgrd
315         WRITE(numout,*)'      oce co2 flux                    = ', TRIM(sn_snd_co2%cldes   ), ' (', TRIM(sn_snd_co2%clcat   ), ')'
316         WRITE(numout,*)'      water level                     = ', TRIM(sn_snd_wlev%cldes  ), ' (', TRIM(sn_snd_wlev%clcat  ), ')' 
317         WRITE(numout,*)'      mean sea level pressure         = ', TRIM(sn_rcv_mslp%cldes  ), ' (', TRIM(sn_rcv_mslp%clcat  ), ')' 
318         WRITE(numout,*)'      surface current to waves        = ', TRIM(sn_snd_crtw%cldes  ), ' (', TRIM(sn_snd_crtw%clcat  ), ')' 
319         WRITE(numout,*)'                      - referential   = ', sn_snd_crtw%clvref 
320         WRITE(numout,*)'                      - orientation   = ', sn_snd_crtw%clvor 
321         WRITE(numout,*)'                      - mesh          = ', sn_snd_crtw%clvgrd 
322         WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel
323         WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask
324      ENDIF
325
326      !                                   ! allocate sbccpl arrays
327      IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )
328     
329      ! ================================ !
330      !   Define the receive interface   !
331      ! ================================ !
332      nrcvinfo(:) = OASIS_idle   ! needed by nrcvinfo(jpr_otx1) if we do not receive ocean stress
333
334      ! for each field: define the OASIS name                              (srcv(:)%clname)
335      !                 define receive or not from the namelist parameters (srcv(:)%laction)
336      !                 define the north fold type of lbc                  (srcv(:)%nsgn)
337
338      ! default definitions of srcv
339      srcv(:)%laction = .FALSE.   ;   srcv(:)%clgrid = 'T'   ;   srcv(:)%nsgn = 1.   ;   srcv(:)%nct = 1
340
341      !                                                      ! ------------------------- !
342      !                                                      ! ice and ocean wind stress !   
343      !                                                      ! ------------------------- !
344      !                                                           ! Name
345      srcv(jpr_otx1)%clname = 'O_OTaux1'      ! 1st ocean component on grid ONE (T or U)
346      srcv(jpr_oty1)%clname = 'O_OTauy1'      ! 2nd   -      -         -     -
347      srcv(jpr_otz1)%clname = 'O_OTauz1'      ! 3rd   -      -         -     -
348      srcv(jpr_otx2)%clname = 'O_OTaux2'      ! 1st ocean component on grid TWO (V)
349      srcv(jpr_oty2)%clname = 'O_OTauy2'      ! 2nd   -      -         -     -
350      srcv(jpr_otz2)%clname = 'O_OTauz2'      ! 3rd   -      -         -     -
351      !
352      srcv(jpr_itx1)%clname = 'O_ITaux1'      ! 1st  ice  component on grid ONE (T, F, I or U)
353      srcv(jpr_ity1)%clname = 'O_ITauy1'      ! 2nd   -      -         -     -
354      srcv(jpr_itz1)%clname = 'O_ITauz1'      ! 3rd   -      -         -     -
355      srcv(jpr_itx2)%clname = 'O_ITaux2'      ! 1st  ice  component on grid TWO (V)
356      srcv(jpr_ity2)%clname = 'O_ITauy2'      ! 2nd   -      -         -     -
357      srcv(jpr_itz2)%clname = 'O_ITauz2'      ! 3rd   -      -         -     -
358      !
359      ! Vectors: change of sign at north fold ONLY if on the local grid
360      IF( TRIM( sn_rcv_tau%cldes ) == 'oce only' .OR. TRIM(sn_rcv_tau%cldes ) == 'oce and ice') THEN ! avoid working with the atmospheric fields if they are not coupled
361      IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' )   srcv(jpr_otx1:jpr_itz2)%nsgn = -1.
362     
363      !                                                           ! Set grid and action
364      SELECT CASE( TRIM( sn_rcv_tau%clvgrd ) )      !  'T', 'U,V', 'U,V,I', 'U,V,F', 'T,I', 'T,F', or 'T,U,V'
365      CASE( 'T' ) 
366         srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
367         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
368         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
369      CASE( 'U,V' ) 
370         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
371         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
372         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'U'        ! ice components given at U-point
373         srcv(jpr_itx2:jpr_itz2)%clgrid  = 'V'        !           and           V-point
374         srcv(jpr_otx1:jpr_itz2)%laction = .TRUE.     ! receive oce and ice components on both grid 1 & 2
375      CASE( 'U,V,T' )
376         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
377         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
378         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'T'        ! ice components given at T-point
379         srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
380         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
381      CASE( 'U,V,I' )
382         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
383         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
384         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'I'        ! ice components given at I-point
385         srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
386         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
387      CASE( 'U,V,F' )
388         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
389         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
390         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'F'        ! ice components given at F-point
391         srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
392         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
393      CASE( 'T,I' ) 
394         srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
395         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'I'        ! ice components given at I-point
396         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
397         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
398      CASE( 'T,F' ) 
399         srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
400         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'F'        ! ice components given at F-point
401         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
402         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
403      CASE( 'T,U,V' )
404         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'T'        ! oce components given at T-point
405         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'U'        ! ice components given at U-point
406         srcv(jpr_itx2:jpr_itz2)%clgrid  = 'V'        !           and           V-point
407         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1 only
408         srcv(jpr_itx1:jpr_itz2)%laction = .TRUE.     ! receive ice components on grid 1 & 2
409      CASE default   
410         CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_tau%clvgrd' )
411      END SELECT
412      !
413      IF( TRIM( sn_rcv_tau%clvref ) == 'spherical' )   &           ! spherical: 3rd component not received
414         &     srcv( (/jpr_otz1, jpr_otz2, jpr_itz1, jpr_itz2/) )%laction = .FALSE. 
415      !
416      IF( TRIM( sn_rcv_tau%clvor  ) == 'local grid' ) THEN        ! already on local grid -> no need of the second grid
417            srcv(jpr_otx2:jpr_otz2)%laction = .FALSE. 
418            srcv(jpr_itx2:jpr_itz2)%laction = .FALSE. 
419            srcv(jpr_oty1)%clgrid = srcv(jpr_oty2)%clgrid   ! not needed but cleaner...
420            srcv(jpr_ity1)%clgrid = srcv(jpr_ity2)%clgrid   ! not needed but cleaner...
421      ENDIF
422      !
423      IF( TRIM( sn_rcv_tau%cldes ) /= 'oce and ice' ) THEN        ! 'oce and ice' case ocean stress on ocean mesh used
424         srcv(jpr_itx1:jpr_itz2)%laction = .FALSE.    ! ice components not received
425         srcv(jpr_itx1)%clgrid = 'U'                  ! ocean stress used after its transformation
426         srcv(jpr_ity1)%clgrid = 'V'                  ! i.e. it is always at U- & V-points for i- & j-comp. resp.
427      ENDIF
428      ENDIF
429
430      !                                                      ! ------------------------- !
431      !                                                      !    freshwater budget      !   E-P
432      !                                                      ! ------------------------- !
433      ! we suppose that atmosphere modele do not make the difference between precipiration (liquide or solid)
434      ! over ice of free ocean within the same atmospheric cell.cd
435      srcv(jpr_rain)%clname = 'OTotRain'      ! Rain = liquid precipitation
436      srcv(jpr_snow)%clname = 'OTotSnow'      ! Snow = solid precipitation
437      srcv(jpr_tevp)%clname = 'OTotEvap'      ! total evaporation (over oce + ice sublimation)
438      srcv(jpr_ievp)%clname = 'OIceEvap'      ! evaporation over ice = sublimation
439      srcv(jpr_sbpr)%clname = 'OSubMPre'      ! sublimation - liquid precipitation - solid precipitation
440      srcv(jpr_semp)%clname = 'OISubMSn'      ! ice solid water budget = sublimation - solid precipitation
441      srcv(jpr_oemp)%clname = 'OOEvaMPr'      ! ocean water budget = ocean Evap - ocean precip
442      SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
443      CASE( 'none'          )       ! nothing to do
444      CASE( 'oce only'      )   ;   srcv(                                 jpr_oemp   )%laction = .TRUE. 
445      CASE( 'conservative'  )
446         srcv( (/jpr_rain, jpr_snow, jpr_ievp, jpr_tevp/) )%laction = .TRUE.
447         IF ( k_ice <= 1 )  srcv(jpr_ievp)%laction = .FALSE.
448      CASE( 'oce and ice'   )   ;   srcv( (/jpr_ievp, jpr_sbpr, jpr_semp, jpr_oemp/) )%laction = .TRUE.
449      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_emp%cldes' )
450      END SELECT
451      !
452      !                                                      ! ------------------------- !
453      !                                                      !     Runoffs & Calving     !   
454      !                                                      ! ------------------------- !
455      srcv(jpr_rnf   )%clname = 'O_Runoff'
456      IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' ) THEN
457         srcv(jpr_rnf)%laction = .TRUE.
458         l_rnfcpl              = .TRUE.                      ! -> no need to read runoffs in sbcrnf
459         ln_rnf                = nn_components /= jp_iam_sas ! -> force to go through sbcrnf if not sas
460         IF(lwp) WRITE(numout,*)
461         IF(lwp) WRITE(numout,*) '   runoffs received from oasis -> force ln_rnf = ', ln_rnf
462      ENDIF
463      !
464      srcv(jpr_cal)%clname = 'OCalving'   ;  IF( TRIM( sn_rcv_cal%cldes) == 'coupled' )   srcv(jpr_cal)%laction = .TRUE.
465      srcv(jpr_isf)%clname = 'OIcshelf'   ;  IF( TRIM( sn_rcv_isf%cldes) == 'coupled' )   srcv(jpr_isf)%laction = .TRUE.
466      srcv(jpr_icb)%clname = 'OIceberg'   ;  IF( TRIM( sn_rcv_icb%cldes) == 'coupled' )   srcv(jpr_icb)%laction = .TRUE.
467
468      IF( srcv(jpr_isf)%laction .AND. ln_isf ) THEN
469         l_isfcpl             = .TRUE.                      ! -> no need to read isf in sbcisf
470         IF(lwp) WRITE(numout,*)
471         IF(lwp) WRITE(numout,*) '   iceshelf received from oasis '
472      ENDIF
473      !
474      !                                                      ! ------------------------- !
475      !                                                      !    non solar radiation    !   Qns
476      !                                                      ! ------------------------- !
477      srcv(jpr_qnsoce)%clname = 'O_QnsOce'
478      srcv(jpr_qnsice)%clname = 'O_QnsIce'
479      srcv(jpr_qnsmix)%clname = 'O_QnsMix'
480      SELECT CASE( TRIM( sn_rcv_qns%cldes ) )
481      CASE( 'none'          )       ! nothing to do
482      CASE( 'oce only'      )   ;   srcv(               jpr_qnsoce   )%laction = .TRUE.
483      CASE( 'conservative'  )   ;   srcv( (/jpr_qnsice, jpr_qnsmix/) )%laction = .TRUE.
484      CASE( 'oce and ice'   )   ;   srcv( (/jpr_qnsice, jpr_qnsoce/) )%laction = .TRUE.
485      CASE( 'mixed oce-ice' )   ;   srcv(               jpr_qnsmix   )%laction = .TRUE. 
486      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_qns%cldes' )
487      END SELECT
488      IF( TRIM( sn_rcv_qns%cldes ) == 'mixed oce-ice' .AND. jpl > 1 ) &
489         CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qns%cldes not currently allowed to be mixed oce-ice for multi-category ice' )
490      !                                                      ! ------------------------- !
491      !                                                      !    solar radiation        !   Qsr
492      !                                                      ! ------------------------- !
493      srcv(jpr_qsroce)%clname = 'O_QsrOce'
494      srcv(jpr_qsrice)%clname = 'O_QsrIce'
495      srcv(jpr_qsrmix)%clname = 'O_QsrMix'
496      SELECT CASE( TRIM( sn_rcv_qsr%cldes ) )
497      CASE( 'none'          )       ! nothing to do
498      CASE( 'oce only'      )   ;   srcv(               jpr_qsroce   )%laction = .TRUE.
499      CASE( 'conservative'  )   ;   srcv( (/jpr_qsrice, jpr_qsrmix/) )%laction = .TRUE.
500      CASE( 'oce and ice'   )   ;   srcv( (/jpr_qsrice, jpr_qsroce/) )%laction = .TRUE.
501      CASE( 'mixed oce-ice' )   ;   srcv(               jpr_qsrmix   )%laction = .TRUE. 
502      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_qsr%cldes' )
503      END SELECT
504      IF( TRIM( sn_rcv_qsr%cldes ) == 'mixed oce-ice' .AND. jpl > 1 ) &
505         CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qsr%cldes not currently allowed to be mixed oce-ice for multi-category ice' )
506      !                                                      ! ------------------------- !
507      !                                                      !   non solar sensitivity   !   d(Qns)/d(T)
508      !                                                      ! ------------------------- !
509      srcv(jpr_dqnsdt)%clname = 'O_dQnsdT'   
510      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'coupled' )   srcv(jpr_dqnsdt)%laction = .TRUE.
511      !
512      ! non solar sensitivity mandatory for LIM ice model
513      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'none' .AND. k_ice /= 0 .AND. k_ice /= 4 .AND. nn_components /= jp_iam_sas ) &
514         CALL ctl_stop( 'sbc_cpl_init: sn_rcv_dqnsdt%cldes must be coupled in namsbc_cpl namelist' )
515      ! non solar sensitivity mandatory for mixed oce-ice solar radiation coupling technique
516      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'none' .AND. TRIM( sn_rcv_qns%cldes ) == 'mixed oce-ice' ) &
517         CALL ctl_stop( 'sbc_cpl_init: namsbc_cpl namelist mismatch between sn_rcv_qns%cldes and sn_rcv_dqnsdt%cldes' )
518      !                                                      ! ------------------------- !
519      !                                                      !      10m wind module      !   
520      !                                                      ! ------------------------- !
521      srcv(jpr_w10m)%clname = 'O_Wind10'   ;   IF( TRIM(sn_rcv_w10m%cldes  ) == 'coupled' )   srcv(jpr_w10m)%laction = .TRUE. 
522      !
523      !                                                      ! ------------------------- !
524      !                                                      !   wind stress module      !   
525      !                                                      ! ------------------------- !
526      srcv(jpr_taum)%clname = 'O_TauMod'   ;   IF( TRIM(sn_rcv_taumod%cldes) == 'coupled' )   srcv(jpr_taum)%laction = .TRUE.
527      lhftau = srcv(jpr_taum)%laction
528
529      !                                                      ! ------------------------- !
530      !                                                      !      Atmospheric CO2      !
531      !                                                      ! ------------------------- !
532      srcv(jpr_co2 )%clname = 'O_AtmCO2'   
533      IF( TRIM(sn_rcv_co2%cldes   ) == 'coupled' )  THEN
534         srcv(jpr_co2 )%laction = .TRUE.
535         l_co2cpl = .TRUE.
536         IF(lwp) WRITE(numout,*)
537         IF(lwp) WRITE(numout,*) '   Atmospheric pco2 received from oasis '
538         IF(lwp) WRITE(numout,*)
539      ENDIF
540
541      !                                                      ! ------------------------- !
542      !                                                      ! Mean Sea Level Pressure   !
543      !                                                      ! ------------------------- !
544      srcv(jpr_mslp)%clname = 'O_MSLP'     ;   IF( TRIM(sn_rcv_mslp%cldes  ) == 'coupled' )    srcv(jpr_mslp)%laction = .TRUE. 
545
546      !                                                      ! ------------------------- !
547      !                                                      !   topmelt and botmelt     !   
548      !                                                      ! ------------------------- !
549      srcv(jpr_topm )%clname = 'OTopMlt'
550      srcv(jpr_botm )%clname = 'OBotMlt'
551      IF( TRIM(sn_rcv_iceflx%cldes) == 'coupled' ) THEN
552         IF ( TRIM( sn_rcv_iceflx%clcat ) == 'yes' ) THEN
553            srcv(jpr_topm:jpr_botm)%nct = jpl
554         ELSE
555            CALL ctl_stop( 'sbc_cpl_init: sn_rcv_iceflx%clcat should always be set to yes currently' )
556         ENDIF
557         srcv(jpr_topm:jpr_botm)%laction = .TRUE.
558      ENDIF
559      !                                                      ! ------------------------- !
560      !                                                      !      Wave breaking        !   
561      !                                                      ! ------------------------- !
562      srcv(jpr_hsig)%clname  = 'O_Hsigwa'    ! significant wave height
563      IF( TRIM(sn_rcv_hsig%cldes  ) == 'coupled' )  THEN
564         srcv(jpr_hsig)%laction = .TRUE.
565         cpl_hsig = .TRUE.
566      ENDIF
567      srcv(jpr_phioc)%clname = 'O_PhiOce'    ! wave to ocean energy
568      IF( TRIM(sn_rcv_phioc%cldes ) == 'coupled' )  THEN
569         srcv(jpr_phioc)%laction = .TRUE.
570         cpl_phioc = .TRUE.
571      ENDIF
572      srcv(jpr_sdrftx)%clname = 'O_Sdrfx'    ! Stokes drift in the u direction
573      IF( TRIM(sn_rcv_sdrfx%cldes ) == 'coupled' )  THEN
574         srcv(jpr_sdrftx)%laction = .TRUE.
575         cpl_sdrftx = .TRUE.
576      ENDIF
577      srcv(jpr_sdrfty)%clname = 'O_Sdrfy'    ! Stokes drift in the v direction
578      IF( TRIM(sn_rcv_sdrfy%cldes ) == 'coupled' )  THEN
579         srcv(jpr_sdrfty)%laction = .TRUE.
580         cpl_sdrfty = .TRUE.
581      ENDIF
582      srcv(jpr_wper)%clname = 'O_WPer'       ! mean wave period
583      IF( TRIM(sn_rcv_wper%cldes  ) == 'coupled' )  THEN
584         srcv(jpr_wper)%laction = .TRUE.
585         cpl_wper = .TRUE.
586      ENDIF
587      srcv(jpr_wfreq)%clname = 'O_WFreq'     ! wave peak frequency
588      IF( TRIM(sn_rcv_wfreq%cldes ) == 'coupled' )  THEN
589         srcv(jpr_wfreq)%laction = .TRUE.
590         cpl_wfreq = .TRUE.
591      ENDIF
592      srcv(jpr_wnum)%clname = 'O_WNum'       ! mean wave number
593      IF( TRIM(sn_rcv_wnum%cldes ) == 'coupled' )  THEN
594         srcv(jpr_wnum)%laction = .TRUE.
595         cpl_wnum = .TRUE.
596      ENDIF
597      srcv(jpr_tauoc)%clname = 'O_TauOce'    ! stress fraction adsorbed by the wave
598      IF( TRIM(sn_rcv_tauoc%cldes ) == 'coupled' )  THEN
599         srcv(jpr_tauoc)%laction = .TRUE.
600         cpl_tauoc = .TRUE.
601      ENDIF
602      srcv(jpr_tauwx)%clname = 'O_Tauwx'      ! ocean stress from wave in the x direction
603      srcv(jpr_tauwy)%clname = 'O_Tauwy'      ! ocean stress from wave in the y direction
604      IF( TRIM(sn_rcv_tauw%cldes ) == 'coupled' )  THEN
605         srcv(jpr_tauwx)%laction = .TRUE.
606         srcv(jpr_tauwy)%laction = .TRUE.
607         cpl_tauw = .TRUE.
608      ENDIF
609      srcv(jpr_wdrag)%clname = 'O_WDrag'     ! neutral surface drag coefficient
610      IF( TRIM(sn_rcv_wdrag%cldes ) == 'coupled' )  THEN
611         srcv(jpr_wdrag)%laction = .TRUE.
612         cpl_wdrag = .TRUE.
613      ENDIF
614      !
615      IF( srcv(jpr_tauoc)%laction .AND. srcv(jpr_tauwx)%laction .AND. srcv(jpr_tauwy)%laction ) &
616            CALL ctl_stop( 'More than one method for modifying the ocean stress has been selected ', &
617                                     '(sn_rcv_tauoc=coupled and sn_rcv_tauw=coupled)' )
618      !
619      !                                                      ! ------------------------------- !
620      !                                                      !   OPA-SAS coupling - rcv by opa !   
621      !                                                      ! ------------------------------- !
622      srcv(jpr_sflx)%clname = 'O_SFLX'
623      srcv(jpr_fice)%clname = 'RIceFrc'
624      !
625      IF( nn_components == jp_iam_opa ) THEN    ! OPA coupled to SAS via OASIS: force received field by OPA (sent by SAS)
626         srcv(:)%laction = .FALSE.   ! force default definition in case of opa <-> sas coupling
627         srcv(:)%clgrid  = 'T'       ! force default definition in case of opa <-> sas coupling
628         srcv(:)%nsgn    = 1.        ! force default definition in case of opa <-> sas coupling
629         srcv( (/jpr_qsroce, jpr_qnsoce, jpr_oemp, jpr_sflx, jpr_fice, jpr_otx1, jpr_oty1, jpr_taum/) )%laction = .TRUE.
630         srcv(jpr_otx1)%clgrid = 'U'        ! oce components given at U-point
631         srcv(jpr_oty1)%clgrid = 'V'        !           and           V-point
632         ! Vectors: change of sign at north fold ONLY if on the local grid
633         srcv( (/jpr_otx1,jpr_oty1/) )%nsgn = -1.
634         sn_rcv_tau%clvgrd = 'U,V'
635         sn_rcv_tau%clvor = 'local grid'
636         sn_rcv_tau%clvref = 'spherical'
637         sn_rcv_emp%cldes = 'oce only'
638         !
639         IF(lwp) THEN                        ! control print
640            WRITE(numout,*)
641            WRITE(numout,*)'               Special conditions for SAS-OPA coupling  '
642            WRITE(numout,*)'               OPA component  '
643            WRITE(numout,*)
644            WRITE(numout,*)'  received fields from SAS component '
645            WRITE(numout,*)'                  ice cover '
646            WRITE(numout,*)'                  oce only EMP  '
647            WRITE(numout,*)'                  salt flux  '
648            WRITE(numout,*)'                  mixed oce-ice solar flux  '
649            WRITE(numout,*)'                  mixed oce-ice non solar flux  '
650            WRITE(numout,*)'                  wind stress U,V on local grid and sperical coordinates '
651            WRITE(numout,*)'                  wind stress module'
652            WRITE(numout,*)
653         ENDIF
654      ENDIF
655      !                                                      ! -------------------------------- !
656      !                                                      !   OPA-SAS coupling - rcv by sas  !   
657      !                                                      ! -------------------------------- !
658      srcv(jpr_toce  )%clname = 'I_SSTSST'
659      srcv(jpr_soce  )%clname = 'I_SSSal'
660      srcv(jpr_ocx1  )%clname = 'I_OCurx1'
661      srcv(jpr_ocy1  )%clname = 'I_OCury1'
662      srcv(jpr_ssh   )%clname = 'I_SSHght'
663      srcv(jpr_e3t1st)%clname = 'I_E3T1st'   
664      srcv(jpr_fraqsr)%clname = 'I_FraQsr'   
665      !
666      IF( nn_components == jp_iam_sas ) THEN
667         IF( .NOT. ln_cpl ) srcv(:)%laction = .FALSE.   ! force default definition in case of opa <-> sas coupling
668         IF( .NOT. ln_cpl ) srcv(:)%clgrid  = 'T'       ! force default definition in case of opa <-> sas coupling
669         IF( .NOT. ln_cpl ) srcv(:)%nsgn    = 1.        ! force default definition in case of opa <-> sas coupling
670         srcv( (/jpr_toce, jpr_soce, jpr_ssh, jpr_fraqsr, jpr_ocx1, jpr_ocy1/) )%laction = .TRUE.
671         srcv( jpr_e3t1st )%laction = .NOT.ln_linssh
672         srcv(jpr_ocx1)%clgrid = 'U'        ! oce components given at U-point
673         srcv(jpr_ocy1)%clgrid = 'V'        !           and           V-point
674         ! Vectors: change of sign at north fold ONLY if on the local grid
675         srcv(jpr_ocx1:jpr_ocy1)%nsgn = -1.
676         ! Change first letter to couple with atmosphere if already coupled OPA
677         ! this is nedeed as each variable name used in the namcouple must be unique:
678         ! for example O_Runoff received by OPA from SAS and therefore O_Runoff received by SAS from the Atmosphere
679         DO jn = 1, jprcv
680            IF ( srcv(jn)%clname(1:1) == "O" ) srcv(jn)%clname = "S"//srcv(jn)%clname(2:LEN(srcv(jn)%clname))
681         END DO
682         !
683         IF(lwp) THEN                        ! control print
684            WRITE(numout,*)
685            WRITE(numout,*)'               Special conditions for SAS-OPA coupling  '
686            WRITE(numout,*)'               SAS component  '
687            WRITE(numout,*)
688            IF( .NOT. ln_cpl ) THEN
689               WRITE(numout,*)'  received fields from OPA component '
690            ELSE
691               WRITE(numout,*)'  Additional received fields from OPA component : '
692            ENDIF
693            WRITE(numout,*)'               sea surface temperature (Celsius) '
694            WRITE(numout,*)'               sea surface salinity ' 
695            WRITE(numout,*)'               surface currents ' 
696            WRITE(numout,*)'               sea surface height ' 
697            WRITE(numout,*)'               thickness of first ocean T level '       
698            WRITE(numout,*)'               fraction of solar net radiation absorbed in the first ocean level'
699            WRITE(numout,*)
700         ENDIF
701      ENDIF
702     
703      ! =================================================== !
704      ! Allocate all parts of frcv used for received fields !
705      ! =================================================== !
706      DO jn = 1, jprcv
707         IF ( srcv(jn)%laction ) ALLOCATE( frcv(jn)%z3(jpi,jpj,srcv(jn)%nct) )
708      END DO
709      ! Allocate taum part of frcv which is used even when not received as coupling field
710      IF ( .NOT. srcv(jpr_taum)%laction ) ALLOCATE( frcv(jpr_taum)%z3(jpi,jpj,srcv(jpr_taum)%nct) )
711      ! Allocate w10m part of frcv which is used even when not received as coupling field
712      IF ( .NOT. srcv(jpr_w10m)%laction ) ALLOCATE( frcv(jpr_w10m)%z3(jpi,jpj,srcv(jpr_w10m)%nct) )
713      ! Allocate jpr_otx1 part of frcv which is used even when not received as coupling field
714      IF ( .NOT. srcv(jpr_otx1)%laction ) ALLOCATE( frcv(jpr_otx1)%z3(jpi,jpj,srcv(jpr_otx1)%nct) )
715      IF ( .NOT. srcv(jpr_oty1)%laction ) ALLOCATE( frcv(jpr_oty1)%z3(jpi,jpj,srcv(jpr_oty1)%nct) )
716      ! Allocate itx1 and ity1 as they are used in sbc_cpl_ice_tau even if srcv(jpr_itx1)%laction = .FALSE.
717      IF( k_ice /= 0 ) THEN
718         IF ( .NOT. srcv(jpr_itx1)%laction ) ALLOCATE( frcv(jpr_itx1)%z3(jpi,jpj,srcv(jpr_itx1)%nct) )
719         IF ( .NOT. srcv(jpr_ity1)%laction ) ALLOCATE( frcv(jpr_ity1)%z3(jpi,jpj,srcv(jpr_ity1)%nct) )
720      END IF
721
722      ! ================================ !
723      !     Define the send interface    !
724      ! ================================ !
725      ! for each field: define the OASIS name                           (ssnd(:)%clname)
726      !                 define send or not from the namelist parameters (ssnd(:)%laction)
727      !                 define the north fold type of lbc               (ssnd(:)%nsgn)
728     
729      ! default definitions of nsnd
730      ssnd(:)%laction = .FALSE.   ;   ssnd(:)%clgrid = 'T'   ;   ssnd(:)%nsgn = 1.  ; ssnd(:)%nct = 1
731         
732      !                                                      ! ------------------------- !
733      !                                                      !    Surface temperature    !
734      !                                                      ! ------------------------- !
735      ssnd(jps_toce)%clname = 'O_SSTSST'
736      ssnd(jps_tice)%clname = 'O_TepIce'
737      ssnd(jps_tmix)%clname = 'O_TepMix'
738      SELECT CASE( TRIM( sn_snd_temp%cldes ) )
739      CASE( 'none'                                 )       ! nothing to do
740      CASE( 'oce only'                             )   ;   ssnd( jps_toce )%laction = .TRUE.
741      CASE( 'oce and ice' , 'weighted oce and ice' )
742         ssnd( (/jps_toce, jps_tice/) )%laction = .TRUE.
743         IF ( TRIM( sn_snd_temp%clcat ) == 'yes' )  ssnd(jps_tice)%nct = jpl
744      CASE( 'mixed oce-ice'                        )   ;   ssnd( jps_tmix )%laction = .TRUE.
745      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_temp%cldes' )
746      END SELECT
747           
748      !                                                      ! ------------------------- !
749      !                                                      !          Albedo           !
750      !                                                      ! ------------------------- !
751      ssnd(jps_albice)%clname = 'O_AlbIce' 
752      ssnd(jps_albmix)%clname = 'O_AlbMix'
753      SELECT CASE( TRIM( sn_snd_alb%cldes ) )
754      CASE( 'none'                 )     ! nothing to do
755      CASE( 'ice' , 'weighted ice' )   ; ssnd(jps_albice)%laction = .TRUE.
756      CASE( 'mixed oce-ice'        )   ; ssnd(jps_albmix)%laction = .TRUE.
757      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_alb%cldes' )
758      END SELECT
759      !
760      ! Need to calculate oceanic albedo if
761      !     1. sending mixed oce-ice albedo or
762      !     2. receiving mixed oce-ice solar radiation
763      IF ( TRIM ( sn_snd_alb%cldes ) == 'mixed oce-ice' .OR. TRIM ( sn_rcv_qsr%cldes ) == 'mixed oce-ice' ) THEN
764         CALL albedo_oce( zaos, zacs )
765         ! Due to lack of information on nebulosity : mean clear/overcast sky
766         albedo_oce_mix(:,:) = ( zacs(:,:) + zaos(:,:) ) * 0.5
767      ENDIF
768
769      !                                                      ! ------------------------- !
770      !                                                      !  Ice fraction & Thickness !
771      !                                                      ! ------------------------- !
772      ssnd(jps_fice)%clname = 'OIceFrc'
773      ssnd(jps_ficet)%clname = 'OIceFrcT' 
774      ssnd(jps_hice)%clname = 'OIceTck'
775      ssnd(jps_hsnw)%clname = 'OSnwTck'
776      IF( k_ice /= 0 ) THEN
777         ssnd(jps_fice)%laction = .TRUE.                  ! if ice treated in the ocean (even in climato case)
778! Currently no namelist entry to determine sending of multi-category ice fraction so use the thickness entry for now
779         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_fice)%nct = jpl
780      ENDIF
781     
782      IF (TRIM( sn_snd_ifrac%cldes )  == 'coupled') ssnd(jps_ficet)%laction = .TRUE. 
783
784      SELECT CASE ( TRIM( sn_snd_thick%cldes ) )
785      CASE( 'none'         )       ! nothing to do
786      CASE( 'ice and snow' ) 
787         ssnd(jps_hice:jps_hsnw)%laction = .TRUE.
788         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) THEN
789            ssnd(jps_hice:jps_hsnw)%nct = jpl
790         ENDIF
791      CASE ( 'weighted ice and snow' ) 
792         ssnd(jps_hice:jps_hsnw)%laction = .TRUE.
793         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_hice:jps_hsnw)%nct = jpl
794      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_thick%cldes' )
795      END SELECT
796
797      !                                                      ! ------------------------- !
798      !                                                      !      Surface current      !
799      !                                                      ! ------------------------- !
800      !        ocean currents              !            ice velocities
801      ssnd(jps_ocx1)%clname = 'O_OCurx1'   ;   ssnd(jps_ivx1)%clname = 'O_IVelx1'
802      ssnd(jps_ocy1)%clname = 'O_OCury1'   ;   ssnd(jps_ivy1)%clname = 'O_IVely1'
803      ssnd(jps_ocz1)%clname = 'O_OCurz1'   ;   ssnd(jps_ivz1)%clname = 'O_IVelz1'
804      ssnd(jps_ocxw)%clname = 'O_OCurxw' 
805      ssnd(jps_ocyw)%clname = 'O_OCuryw' 
806      !
807      ssnd(jps_ocx1:jps_ivz1)%nsgn = -1.   ! vectors: change of the sign at the north fold
808
809      IF( sn_snd_crt%clvgrd == 'U,V' ) THEN
810         ssnd(jps_ocx1)%clgrid = 'U' ; ssnd(jps_ocy1)%clgrid = 'V'
811      ELSE IF( sn_snd_crt%clvgrd /= 'T' ) THEN 
812         CALL ctl_stop( 'sn_snd_crt%clvgrd must be equal to T' )
813         ssnd(jps_ocx1:jps_ivz1)%clgrid  = 'T'      ! all oce and ice components on the same unique grid
814      ENDIF
815      ssnd(jps_ocx1:jps_ivz1)%laction = .TRUE.   ! default: all are send
816      IF( TRIM( sn_snd_crt%clvref ) == 'spherical' )   ssnd( (/jps_ocz1, jps_ivz1/) )%laction = .FALSE. 
817      IF( TRIM( sn_snd_crt%clvor ) == 'eastward-northward' ) ssnd(jps_ocx1:jps_ivz1)%nsgn = 1.
818      SELECT CASE( TRIM( sn_snd_crt%cldes ) )
819      CASE( 'none'                 )   ;   ssnd(jps_ocx1:jps_ivz1)%laction = .FALSE.
820      CASE( 'oce only'             )   ;   ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE.
821      CASE( 'weighted oce and ice' )   !   nothing to do
822      CASE( 'mixed oce-ice'        )   ;   ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE.
823      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_crt%cldes' )
824      END SELECT
825
826      ssnd(jps_ocxw:jps_ocyw)%nsgn = -1.   ! vectors: change of the sign at the north fold
827       
828      IF( sn_snd_crtw%clvgrd == 'U,V' ) THEN
829         ssnd(jps_ocxw)%clgrid = 'U' ; ssnd(jps_ocyw)%clgrid = 'V' 
830      ELSE IF( sn_snd_crtw%clvgrd /= 'T' ) THEN
831         CALL ctl_stop( 'sn_snd_crtw%clvgrd must be equal to T' ) 
832      ENDIF
833      IF( TRIM( sn_snd_crtw%clvor ) == 'eastward-northward' ) ssnd(jps_ocxw:jps_ocyw)%nsgn = 1. 
834      SELECT CASE( TRIM( sn_snd_crtw%cldes ) ) 
835         CASE( 'none'                 )   ; ssnd(jps_ocxw:jps_ocyw)%laction = .FALSE. 
836         CASE( 'oce only'             )   ; ssnd(jps_ocxw:jps_ocyw)%laction = .TRUE. 
837         CASE( 'weighted oce and ice' )   !   nothing to do
838         CASE( 'mixed oce-ice'        )   ; ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE. 
839         CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_crtw%cldes' ) 
840      END SELECT 
841
842      !                                                      ! ------------------------- !
843      !                                                      !          CO2 flux         !
844      !                                                      ! ------------------------- !
845      ssnd(jps_co2)%clname = 'O_CO2FLX' ;  IF( TRIM(sn_snd_co2%cldes) == 'coupled' )    ssnd(jps_co2 )%laction = .TRUE.
846
847      !                                                      ! ------------------------- !
848      !                                                      !     Sea surface height    !
849      !                                                      ! ------------------------- !
850      ssnd(jps_wlev)%clname = 'O_Wlevel' ;  IF( TRIM(sn_snd_wlev%cldes) == 'coupled' )   ssnd(jps_wlev)%laction = .TRUE. 
851
852      !                                                      ! ------------------------------- !
853      !                                                      !   OPA-SAS coupling - snd by opa !   
854      !                                                      ! ------------------------------- !
855      ssnd(jps_ssh   )%clname = 'O_SSHght' 
856      ssnd(jps_soce  )%clname = 'O_SSSal' 
857      ssnd(jps_e3t1st)%clname = 'O_E3T1st'   
858      ssnd(jps_fraqsr)%clname = 'O_FraQsr'
859      !
860      IF( nn_components == jp_iam_opa ) THEN
861         ssnd(:)%laction = .FALSE.   ! force default definition in case of opa <-> sas coupling
862         ssnd( (/jps_toce, jps_soce, jps_ssh, jps_fraqsr, jps_ocx1, jps_ocy1/) )%laction = .TRUE.
863         ssnd( jps_e3t1st )%laction = .NOT.ln_linssh
864         ! vector definition: not used but cleaner...
865         ssnd(jps_ocx1)%clgrid  = 'U'        ! oce components given at U-point
866         ssnd(jps_ocy1)%clgrid  = 'V'        !           and           V-point
867         sn_snd_crt%clvgrd = 'U,V'
868         sn_snd_crt%clvor = 'local grid'
869         sn_snd_crt%clvref = 'spherical'
870         !
871         IF(lwp) THEN                        ! control print
872            WRITE(numout,*)
873            WRITE(numout,*)'  sent fields to SAS component '
874            WRITE(numout,*)'               sea surface temperature (T before, Celsius) '
875            WRITE(numout,*)'               sea surface salinity ' 
876            WRITE(numout,*)'               surface currents U,V on local grid and spherical coordinates' 
877            WRITE(numout,*)'               sea surface height ' 
878            WRITE(numout,*)'               thickness of first ocean T level '       
879            WRITE(numout,*)'               fraction of solar net radiation absorbed in the first ocean level'
880            WRITE(numout,*)
881         ENDIF
882      ENDIF
883      !                                                      ! ------------------------------- !
884      !                                                      !   OPA-SAS coupling - snd by sas !   
885      !                                                      ! ------------------------------- !
886      ssnd(jps_sflx  )%clname = 'I_SFLX'     
887      ssnd(jps_fice2 )%clname = 'IIceFrc'
888      ssnd(jps_qsroce)%clname = 'I_QsrOce'   
889      ssnd(jps_qnsoce)%clname = 'I_QnsOce'   
890      ssnd(jps_oemp  )%clname = 'IOEvaMPr' 
891      ssnd(jps_otx1  )%clname = 'I_OTaux1'   
892      ssnd(jps_oty1  )%clname = 'I_OTauy1'   
893      ssnd(jps_rnf   )%clname = 'I_Runoff'   
894      ssnd(jps_taum  )%clname = 'I_TauMod'   
895      !
896      IF( nn_components == jp_iam_sas ) THEN
897         IF( .NOT. ln_cpl ) ssnd(:)%laction = .FALSE.   ! force default definition in case of opa <-> sas coupling
898         ssnd( (/jps_qsroce, jps_qnsoce, jps_oemp, jps_fice2, jps_sflx, jps_otx1, jps_oty1, jps_taum/) )%laction = .TRUE.
899         !
900         ! Change first letter to couple with atmosphere if already coupled with sea_ice
901         ! this is nedeed as each variable name used in the namcouple must be unique:
902         ! for example O_SSTSST sent by OPA to SAS and therefore S_SSTSST sent by SAS to the Atmosphere
903         DO jn = 1, jpsnd
904            IF ( ssnd(jn)%clname(1:1) == "O" ) ssnd(jn)%clname = "S"//ssnd(jn)%clname(2:LEN(ssnd(jn)%clname))
905         END DO
906         !
907         IF(lwp) THEN                        ! control print
908            WRITE(numout,*)
909            IF( .NOT. ln_cpl ) THEN
910               WRITE(numout,*)'  sent fields to OPA component '
911            ELSE
912               WRITE(numout,*)'  Additional sent fields to OPA component : '
913            ENDIF
914            WRITE(numout,*)'                  ice cover '
915            WRITE(numout,*)'                  oce only EMP  '
916            WRITE(numout,*)'                  salt flux  '
917            WRITE(numout,*)'                  mixed oce-ice solar flux  '
918            WRITE(numout,*)'                  mixed oce-ice non solar flux  '
919            WRITE(numout,*)'                  wind stress U,V components'
920            WRITE(numout,*)'                  wind stress module'
921         ENDIF
922      ENDIF
923
924      !
925      ! ================================ !
926      !   initialisation of the coupler  !
927      ! ================================ !
928
929      CALL cpl_define(jprcv, jpsnd, nn_cplmodel)
930     
931      IF (ln_usecplmask) THEN
932         xcplmask(:,:,:) = 0.
933         CALL iom_open( 'cplmask', inum )
934         CALL iom_get( inum, jpdom_unknown, 'cplmask', xcplmask(1:nlci,1:nlcj,1:nn_cplmodel),   &
935            &          kstart = (/ mig(1),mjg(1),1 /), kcount = (/ nlci,nlcj,nn_cplmodel /) )
936         CALL iom_close( inum )
937      ELSE
938         xcplmask(:,:,:) = 1.
939      ENDIF
940      xcplmask(:,:,0) = 1. - SUM( xcplmask(:,:,1:nn_cplmodel), dim = 3 )
941      !
942      ncpl_qsr_freq = cpl_freq( 'O_QsrOce' ) + cpl_freq( 'O_QsrMix' ) + cpl_freq( 'I_QsrOce' ) + cpl_freq( 'I_QsrMix' )
943      IF( ln_dm2dc .AND. ln_cpl .AND. ncpl_qsr_freq /= 86400 )   &
944         &   CALL ctl_stop( 'sbc_cpl_init: diurnal cycle reconstruction (ln_dm2dc) needs daily couping for solar radiation' )
945      IF( ln_dm2dc .AND. ln_cpl ) ncpl_qsr_freq = 86400 / ncpl_qsr_freq
946
947      CALL wrk_dealloc( jpi,jpj,   zacs, zaos )
948      !
949      IF( nn_timing == 1 )   CALL timing_stop('sbc_cpl_init')
950      !
951   END SUBROUTINE sbc_cpl_init
952
953
954   SUBROUTINE sbc_cpl_rcv( kt, k_fsbc, k_ice )     
955      !!----------------------------------------------------------------------
956      !!             ***  ROUTINE sbc_cpl_rcv  ***
957      !!
958      !! ** Purpose :   provide the stress over the ocean and, if no sea-ice,
959      !!                provide the ocean heat and freshwater fluxes.
960      !!
961      !! ** Method  : - Receive all the atmospheric fields (stored in frcv array). called at each time step.
962      !!                OASIS controls if there is something do receive or not. nrcvinfo contains the info
963      !!                to know if the field was really received or not
964      !!
965      !!              --> If ocean stress was really received:
966      !!
967      !!                  - transform the received ocean stress vector from the received
968      !!                 referential and grid into an atmosphere-ocean stress in
969      !!                 the (i,j) ocean referencial and at the ocean velocity point.
970      !!                    The received stress are :
971      !!                     - defined by 3 components (if cartesian coordinate)
972      !!                            or by 2 components (if spherical)
973      !!                     - oriented along geographical   coordinate (if eastward-northward)
974      !!                            or  along the local grid coordinate (if local grid)
975      !!                     - given at U- and V-point, resp.   if received on 2 grids
976      !!                            or at T-point               if received on 1 grid
977      !!                    Therefore and if necessary, they are successively
978      !!                  processed in order to obtain them
979      !!                     first  as  2 components on the sphere
980      !!                     second as  2 components oriented along the local grid
981      !!                     third  as  2 components on the U,V grid
982      !!
983      !!              -->
984      !!
985      !!              - In 'ocean only' case, non solar and solar ocean heat fluxes
986      !!             and total ocean freshwater fluxes 
987      !!
988      !! ** Method  :   receive all fields from the atmosphere and transform
989      !!              them into ocean surface boundary condition fields
990      !!
991      !! ** Action  :   update  utau, vtau   ocean stress at U,V grid
992      !!                        taum         wind stress module at T-point
993      !!                        wndm         wind speed  module at T-point over free ocean or leads in presence of sea-ice
994      !!                        qns          non solar heat fluxes including emp heat content    (ocean only case)
995      !!                                     and the latent heat flux of solid precip. melting
996      !!                        qsr          solar ocean heat fluxes   (ocean only case)
997      !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case)
998      !!----------------------------------------------------------------------
999      USE zdf_oce,  ONLY : ln_zdfqiao
1000
1001      IMPLICIT NONE
1002
1003      INTEGER, INTENT(in)           ::   kt          ! ocean model time step index
1004      INTEGER, INTENT(in)           ::   k_fsbc      ! frequency of sbc (-> ice model) computation
1005      INTEGER, INTENT(in)           ::   k_ice       ! ice management in the sbc (=0/1/2/3)
1006      !!
1007      LOGICAL  ::   llnewtx, llnewtau      ! update wind stress components and module??
1008      INTEGER  ::   ji, jj, jn             ! dummy loop indices
1009      INTEGER  ::   isec                   ! number of seconds since nit000 (assuming rdt did not change since nit000)
1010      REAL(wp) ::   zcumulneg, zcumulpos   ! temporary scalars     
1011      REAL(wp) ::   zcoef                  ! temporary scalar
1012      REAL(wp) ::   zrhoa  = 1.22          ! Air density kg/m3
1013      REAL(wp) ::   zcdrag = 1.5e-3        ! drag coefficient
1014      REAL(wp) ::   zzx, zzy               ! temporary variables
1015      REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty, zmsk, zemp, zqns, zqsr
1016      !!----------------------------------------------------------------------
1017      !
1018      IF( nn_timing == 1 )   CALL timing_start('sbc_cpl_rcv')
1019      !
1020      CALL wrk_alloc( jpi,jpj,   ztx, zty, zmsk, zemp, zqns, zqsr )
1021      !
1022      IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
1023      !
1024      !                                                      ! ======================================================= !
1025      !                                                      ! Receive all the atmos. fields (including ice information)
1026      !                                                      ! ======================================================= !
1027      isec = ( kt - nit000 ) * NINT( rdt )                      ! date of exchanges
1028      DO jn = 1, jprcv                                          ! received fields sent by the atmosphere
1029         IF( srcv(jn)%laction )   CALL cpl_rcv( jn, isec, frcv(jn)%z3, xcplmask(:,:,1:nn_cplmodel), nrcvinfo(jn) )
1030      END DO
1031
1032      !                                                      ! ========================= !
1033      IF( srcv(jpr_otx1)%laction ) THEN                      !  ocean stress components  !
1034         !                                                   ! ========================= !
1035         ! define frcv(jpr_otx1)%z3(:,:,1) and frcv(jpr_oty1)%z3(:,:,1): stress at U/V point along model grid
1036         ! => need to be done only when we receive the field
1037         IF(  nrcvinfo(jpr_otx1) == OASIS_Rcv ) THEN
1038            !
1039            IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN            ! 2 components on the sphere
1040               !                                                       ! (cartesian to spherical -> 3 to 2 components)
1041               !
1042               CALL geo2oce( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), frcv(jpr_otz1)%z3(:,:,1),   &
1043                  &          srcv(jpr_otx1)%clgrid, ztx, zty )
1044               frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 1st grid
1045               frcv(jpr_oty1)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 1st grid
1046               !
1047               IF( srcv(jpr_otx2)%laction ) THEN
1048                  CALL geo2oce( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), frcv(jpr_otz2)%z3(:,:,1),   &
1049                     &          srcv(jpr_otx2)%clgrid, ztx, zty )
1050                  frcv(jpr_otx2)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 2nd grid
1051                  frcv(jpr_oty2)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 2nd grid
1052               ENDIF
1053               !
1054            ENDIF
1055            !
1056            IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN   ! 2 components oriented along the local grid
1057               !                                                       ! (geographical to local grid -> rotate the components)
1058               CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->i', ztx )   
1059               IF( srcv(jpr_otx2)%laction ) THEN
1060                  CALL rot_rep( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), srcv(jpr_otx2)%clgrid, 'en->j', zty )   
1061               ELSE
1062                  CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->j', zty ) 
1063               ENDIF
1064               frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid
1065               frcv(jpr_oty1)%z3(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 2nd grid
1066            ENDIF
1067            !                             
1068            IF( srcv(jpr_otx1)%clgrid == 'T' ) THEN
1069               DO jj = 2, jpjm1                                          ! T ==> (U,V)
1070                  DO ji = fs_2, fs_jpim1   ! vector opt.
1071                     frcv(jpr_otx1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_otx1)%z3(ji+1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1) )
1072                     frcv(jpr_oty1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_oty1)%z3(ji  ,jj+1,1) + frcv(jpr_oty1)%z3(ji,jj,1) )
1073                  END DO
1074               END DO
1075               CALL lbc_lnk( frcv(jpr_otx1)%z3(:,:,1), 'U',  -1. )   ;   CALL lbc_lnk( frcv(jpr_oty1)%z3(:,:,1), 'V',  -1. )
1076            ENDIF
1077            llnewtx = .TRUE.
1078         ELSE
1079            llnewtx = .FALSE.
1080         ENDIF
1081         !                                                   ! ========================= !
1082      ELSE                                                   !   No dynamical coupling   !
1083         !                                                   ! ========================= !
1084         frcv(jpr_otx1)%z3(:,:,1) = 0.e0                               ! here simply set to zero
1085         frcv(jpr_oty1)%z3(:,:,1) = 0.e0                               ! an external read in a file can be added instead
1086         llnewtx = .TRUE.
1087         !
1088      ENDIF
1089      !                                                      ! ========================= !
1090      !                                                      !    wind stress module     !   (taum)
1091      !                                                      ! ========================= !
1092      !
1093      IF( .NOT. srcv(jpr_taum)%laction ) THEN                    ! compute wind stress module from its components if not received
1094         ! => need to be done only when otx1 was changed
1095         IF( llnewtx ) THEN
1096            DO jj = 2, jpjm1
1097               DO ji = fs_2, fs_jpim1   ! vect. opt.
1098                  zzx = frcv(jpr_otx1)%z3(ji-1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1)
1099                  zzy = frcv(jpr_oty1)%z3(ji  ,jj-1,1) + frcv(jpr_oty1)%z3(ji,jj,1)
1100                  frcv(jpr_taum)%z3(ji,jj,1) = 0.5 * SQRT( zzx * zzx + zzy * zzy )
1101               END DO
1102            END DO
1103            CALL lbc_lnk( frcv(jpr_taum)%z3(:,:,1), 'T', 1. )
1104            llnewtau = .TRUE.
1105         ELSE
1106            llnewtau = .FALSE.
1107         ENDIF
1108      ELSE
1109         llnewtau = nrcvinfo(jpr_taum) == OASIS_Rcv
1110         ! Stress module can be negative when received (interpolation problem)
1111         IF( llnewtau ) THEN
1112            frcv(jpr_taum)%z3(:,:,1) = MAX( 0._wp, frcv(jpr_taum)%z3(:,:,1) )
1113         ENDIF
1114      ENDIF
1115      !
1116      !                                                      ! ========================= !
1117      !                                                      !      10 m wind speed      !   (wndm)
1118      !                                                      ! ========================= !
1119      !
1120      IF( .NOT. srcv(jpr_w10m)%laction ) THEN                    ! compute wind spreed from wind stress module if not received 
1121         ! => need to be done only when taumod was changed
1122         IF( llnewtau ) THEN
1123            zcoef = 1. / ( zrhoa * zcdrag ) 
1124            DO jj = 1, jpj
1125               DO ji = 1, jpi 
1126                  frcv(jpr_w10m)%z3(ji,jj,1) = SQRT( frcv(jpr_taum)%z3(ji,jj,1) * zcoef )
1127               END DO
1128            END DO
1129         ENDIF
1130      ENDIF
1131
1132      ! u(v)tau and taum will be modified by ice model
1133      ! -> need to be reset before each call of the ice/fsbc     
1134      IF( MOD( kt-1, k_fsbc ) == 0 ) THEN
1135         !
1136         IF( ln_mixcpl ) THEN
1137            utau(:,:) = utau(:,:) * xcplmask(:,:,0) + frcv(jpr_otx1)%z3(:,:,1) * zmsk(:,:)
1138            vtau(:,:) = vtau(:,:) * xcplmask(:,:,0) + frcv(jpr_oty1)%z3(:,:,1) * zmsk(:,:)
1139            taum(:,:) = taum(:,:) * xcplmask(:,:,0) + frcv(jpr_taum)%z3(:,:,1) * zmsk(:,:)
1140            wndm(:,:) = wndm(:,:) * xcplmask(:,:,0) + frcv(jpr_w10m)%z3(:,:,1) * zmsk(:,:)
1141         ELSE
1142            utau(:,:) = frcv(jpr_otx1)%z3(:,:,1)
1143            vtau(:,:) = frcv(jpr_oty1)%z3(:,:,1)
1144            taum(:,:) = frcv(jpr_taum)%z3(:,:,1)
1145            wndm(:,:) = frcv(jpr_w10m)%z3(:,:,1)
1146         ENDIF
1147         CALL iom_put( "taum_oce", taum )   ! output wind stress module
1148         
1149      ENDIF
1150
1151      !                                                      ! ================== !
1152      !                                                      ! atmosph. CO2 (ppm) !
1153      !                                                      ! ================== !
1154      IF( srcv(jpr_co2)%laction )   atm_co2(:,:) = frcv(jpr_co2)%z3(:,:,1)
1155      !
1156      !                                                      ! ========================= !
1157      !                                                      ! Mean Sea Level Pressure   !   (taum)
1158      !                                                      ! ========================= !
1159      !
1160      IF( srcv(jpr_mslp)%laction ) THEN                    ! UKMO SHELF effect of atmospheric pressure on SSH
1161          IF( kt /= nit000 )   ssh_ibb(:,:) = ssh_ib(:,:)    !* Swap of ssh_ib fields
1162
1163          r1_grau = 1.e0 / (grav * rau0)               !* constant for optimization
1164          ssh_ib(:,:) = - ( frcv(jpr_mslp)%z3(:,:,1) - rpref ) * r1_grau    ! equivalent ssh (inverse barometer)
1165          apr   (:,:) =     frcv(jpr_mslp)%z3(:,:,1)                         !atmospheric pressure
1166   
1167          IF( kt == nit000 ) ssh_ibb(:,:) = ssh_ib(:,:)  ! correct this later (read from restart if possible)
1168      END IF 
1169      !
1170      IF( ln_sdw ) THEN  ! Stokes Drift correction activated
1171      !                                                      ! ========================= !
1172      !                                                      !       Stokes drift u      !
1173      !                                                      ! ========================= !
1174         IF( srcv(jpr_sdrftx)%laction ) ut0sd(:,:) = frcv(jpr_sdrftx)%z3(:,:,1)
1175      !
1176      !                                                      ! ========================= !
1177      !                                                      !       Stokes drift v      !
1178      !                                                      ! ========================= !
1179         IF( srcv(jpr_sdrfty)%laction ) vt0sd(:,:) = frcv(jpr_sdrfty)%z3(:,:,1)
1180      !
1181      !                                                      ! ========================= !
1182      !                                                      !      Wave mean period     !
1183      !                                                      ! ========================= !
1184         IF( srcv(jpr_wper)%laction ) wmp(:,:) = frcv(jpr_wper)%z3(:,:,1)
1185      !
1186      !                                                      ! ========================= !
1187      !                                                      !  Significant wave height  !
1188      !                                                      ! ========================= !
1189         IF( srcv(jpr_hsig)%laction ) hsw(:,:) = frcv(jpr_hsig)%z3(:,:,1)
1190      !
1191      !                                                      ! ========================= ! 
1192      !                                                      !    Wave peak frequency    !
1193      !                                                      ! ========================= ! 
1194         IF( srcv(jpr_wfreq)%laction ) wfreq(:,:) = frcv(jpr_wfreq)%z3(:,:,1)
1195      !
1196      !                                                      ! ========================= !
1197      !                                                      !    Vertical mixing Qiao   !
1198      !                                                      ! ========================= !
1199         IF( srcv(jpr_wnum)%laction .AND. ln_zdfqiao ) wnum(:,:) = frcv(jpr_wnum)%z3(:,:,1)
1200
1201         ! Calculate the 3D Stokes drift both in coupled and not fully uncoupled mode
1202         IF( srcv(jpr_sdrftx)%laction .OR. srcv(jpr_sdrfty)%laction .OR. srcv(jpr_wper)%laction &
1203                                      .OR. srcv(jpr_hsig)%laction   .OR. srcv(jpr_wfreq)%laction ) THEN
1204            CALL sbc_stokes()
1205         ENDIF
1206      ENDIF
1207      !                                                      ! ========================= !
1208      !                                                      ! Stress adsorbed by waves  !
1209      !                                                      ! ========================= !
1210      IF( srcv(jpr_tauoc)%laction .AND. ln_tauoc ) tauoc_wave(:,:) = frcv(jpr_tauoc)%z3(:,:,1)
1211
1212      !                                                      ! ========================= ! 
1213      !                                                      ! Stress component by waves !
1214      !                                                      ! ========================= ! 
1215      IF( srcv(jpr_tauwx)%laction .AND. srcv(jpr_tauwy)%laction .AND. ln_tauw ) THEN
1216         tauw_x(:,:) = frcv(jpr_tauwx)%z3(:,:,1)
1217         tauw_y(:,:) = frcv(jpr_tauwy)%z3(:,:,1)
1218      ENDIF
1219
1220      !                                                      ! ========================= !
1221      !                                                      !   Wave drag coefficient   !
1222      !                                                      ! ========================= !
1223      IF( srcv(jpr_wdrag)%laction .AND. ln_cdgw ) cdn_wave(:,:) = frcv(jpr_wdrag)%z3(:,:,1)
1224
1225      !  Fields received by SAS when OASIS coupling
1226      !  (arrays no more filled at sbcssm stage)
1227      !                                                      ! ================== !
1228      !                                                      !        SSS         !
1229      !                                                      ! ================== !
1230      IF( srcv(jpr_soce)%laction ) THEN                      ! received by sas in case of opa <-> sas coupling
1231         sss_m(:,:) = frcv(jpr_soce)%z3(:,:,1)
1232         CALL iom_put( 'sss_m', sss_m )
1233      ENDIF
1234      !                                               
1235      !                                                      ! ================== !
1236      !                                                      !        SST         !
1237      !                                                      ! ================== !
1238      IF( srcv(jpr_toce)%laction ) THEN                      ! received by sas in case of opa <-> sas coupling
1239         sst_m(:,:) = frcv(jpr_toce)%z3(:,:,1)
1240         IF( srcv(jpr_soce)%laction .AND. l_useCT ) THEN    ! make sure that sst_m is the potential temperature
1241            sst_m(:,:) = eos_pt_from_ct( sst_m(:,:), sss_m(:,:) )
1242         ENDIF
1243      ENDIF
1244      !                                                      ! ================== !
1245      !                                                      !        SSH         !
1246      !                                                      ! ================== !
1247      IF( srcv(jpr_ssh )%laction ) THEN                      ! received by sas in case of opa <-> sas coupling
1248         ssh_m(:,:) = frcv(jpr_ssh )%z3(:,:,1)
1249         CALL iom_put( 'ssh_m', ssh_m )
1250      ENDIF
1251      !                                                      ! ================== !
1252      !                                                      !  surface currents  !
1253      !                                                      ! ================== !
1254      IF( srcv(jpr_ocx1)%laction ) THEN                      ! received by sas in case of opa <-> sas coupling
1255         ssu_m(:,:) = frcv(jpr_ocx1)%z3(:,:,1)
1256         ub (:,:,1) = ssu_m(:,:)                             ! will be used in sbcice_lim in the call of lim_sbc_tau
1257         un (:,:,1) = ssu_m(:,:)                             ! will be used in sbc_cpl_snd if atmosphere coupling
1258         CALL iom_put( 'ssu_m', ssu_m )
1259      ENDIF
1260      IF( srcv(jpr_ocy1)%laction ) THEN
1261         ssv_m(:,:) = frcv(jpr_ocy1)%z3(:,:,1)
1262         vb (:,:,1) = ssv_m(:,:)                             ! will be used in sbcice_lim in the call of lim_sbc_tau
1263         vn (:,:,1) = ssv_m(:,:)                             ! will be used in sbc_cpl_snd if atmosphere coupling
1264         CALL iom_put( 'ssv_m', ssv_m )
1265      ENDIF
1266      !                                                      ! ======================== !
1267      !                                                      !  first T level thickness !
1268      !                                                      ! ======================== !
1269      IF( srcv(jpr_e3t1st )%laction ) THEN                   ! received by sas in case of opa <-> sas coupling
1270         e3t_m(:,:) = frcv(jpr_e3t1st )%z3(:,:,1)
1271         CALL iom_put( 'e3t_m', e3t_m(:,:) )
1272      ENDIF
1273      !                                                      ! ================================ !
1274      !                                                      !  fraction of solar net radiation !
1275      !                                                      ! ================================ !
1276      IF( srcv(jpr_fraqsr)%laction ) THEN                    ! received by sas in case of opa <-> sas coupling
1277         frq_m(:,:) = frcv(jpr_fraqsr)%z3(:,:,1)
1278         CALL iom_put( 'frq_m', frq_m )
1279      ENDIF
1280     
1281      !                                                      ! ========================= !
1282      IF( k_ice <= 1 .AND. MOD( kt-1, k_fsbc ) == 0 ) THEN   !  heat & freshwater fluxes ! (Ocean only case)
1283         !                                                   ! ========================= !
1284         !
1285         !                                                       ! total freshwater fluxes over the ocean (emp)
1286         IF( srcv(jpr_oemp)%laction .OR. srcv(jpr_rain)%laction ) THEN
1287            SELECT CASE( TRIM( sn_rcv_emp%cldes ) )                                    ! evaporation - precipitation
1288            CASE( 'conservative' )
1289               zemp(:,:) = frcv(jpr_tevp)%z3(:,:,1) - ( frcv(jpr_rain)%z3(:,:,1) + frcv(jpr_snow)%z3(:,:,1) )
1290            CASE( 'oce only', 'oce and ice' )
1291               zemp(:,:) = frcv(jpr_oemp)%z3(:,:,1)
1292            CASE default
1293               CALL ctl_stop( 'sbc_cpl_rcv: wrong definition of sn_rcv_emp%cldes' )
1294            END SELECT
1295         ELSE
1296            zemp(:,:) = 0._wp
1297         ENDIF
1298         !
1299         !                                                        ! runoffs and calving (added in emp)
1300         IF( srcv(jpr_rnf)%laction )     rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
1301         IF( srcv(jpr_cal)%laction )     zemp(:,:) = zemp(:,:) - frcv(jpr_cal)%z3(:,:,1)
1302 
1303         IF( srcv(jpr_icb)%laction )  THEN
1304             fwficb(:,:) = frcv(jpr_icb)%z3(:,:,1)
1305             rnf(:,:)    = rnf(:,:) + fwficb(:,:)   ! iceberg added to runfofs
1306         ENDIF
1307         IF( srcv(jpr_isf)%laction )  fwfisf(:,:) = - frcv(jpr_isf)%z3(:,:,1)  ! fresh water flux from the isf (fwfisf <0 mean melting) 
1308       
1309         IF( ln_mixcpl ) THEN   ;   emp(:,:) = emp(:,:) * xcplmask(:,:,0) + zemp(:,:) * zmsk(:,:)
1310         ELSE                   ;   emp(:,:) =                              zemp(:,:)
1311         ENDIF
1312         !
1313         !                                                       ! non solar heat flux over the ocean (qns)
1314         IF(      srcv(jpr_qnsoce)%laction ) THEN   ;   zqns(:,:) = frcv(jpr_qnsoce)%z3(:,:,1)
1315         ELSE IF( srcv(jpr_qnsmix)%laction ) THEN   ;   zqns(:,:) = frcv(jpr_qnsmix)%z3(:,:,1)
1316         ELSE                                       ;   zqns(:,:) = 0._wp
1317         END IF
1318         ! update qns over the free ocean with:
1319         IF( nn_components /= jp_iam_opa ) THEN
1320            zqns(:,:) =  zqns(:,:) - zemp(:,:) * sst_m(:,:) * rcp         ! remove heat content due to mass flux (assumed to be at SST)
1321            IF( srcv(jpr_snow  )%laction ) THEN
1322               zqns(:,:) = zqns(:,:) - frcv(jpr_snow)%z3(:,:,1) * lfus    ! energy for melting solid precipitation over the free ocean
1323            ENDIF
1324         ENDIF
1325         !
1326         IF( srcv(jpr_icb)%laction )  zqns(:,:) = zqns(:,:) - frcv(jpr_icb)%z3(:,:,1) * lfus ! remove heat content associated to iceberg melting
1327         !
1328         IF( ln_mixcpl ) THEN   ;   qns(:,:) = qns(:,:) * xcplmask(:,:,0) + zqns(:,:) * zmsk(:,:)
1329         ELSE                   ;   qns(:,:) =                              zqns(:,:)
1330         ENDIF
1331
1332         !                                                       ! solar flux over the ocean          (qsr)
1333         IF     ( srcv(jpr_qsroce)%laction ) THEN   ;   zqsr(:,:) = frcv(jpr_qsroce)%z3(:,:,1)
1334         ELSE IF( srcv(jpr_qsrmix)%laction ) then   ;   zqsr(:,:) = frcv(jpr_qsrmix)%z3(:,:,1)
1335         ELSE                                       ;   zqsr(:,:) = 0._wp
1336         ENDIF
1337         IF( ln_dm2dc .AND. ln_cpl )   zqsr(:,:) = sbc_dcy( zqsr )   ! modify qsr to include the diurnal cycle
1338         IF( ln_mixcpl ) THEN   ;   qsr(:,:) = qsr(:,:) * xcplmask(:,:,0) + zqsr(:,:) * zmsk(:,:)
1339         ELSE                   ;   qsr(:,:) =                              zqsr(:,:)
1340         ENDIF
1341         !
1342         ! salt flux over the ocean (received by opa in case of opa <-> sas coupling)
1343         IF( srcv(jpr_sflx )%laction )   sfx(:,:) = frcv(jpr_sflx  )%z3(:,:,1)
1344         ! Ice cover  (received by opa in case of opa <-> sas coupling)
1345         IF( srcv(jpr_fice )%laction )   fr_i(:,:) = frcv(jpr_fice )%z3(:,:,1)
1346         !
1347      ENDIF
1348      !
1349      CALL wrk_dealloc( jpi,jpj,   ztx, zty, zmsk, zemp, zqns, zqsr )
1350      !
1351      IF( nn_timing == 1 )   CALL timing_stop('sbc_cpl_rcv')
1352      !
1353   END SUBROUTINE sbc_cpl_rcv
1354   
1355
1356   SUBROUTINE sbc_cpl_ice_tau( p_taui, p_tauj )     
1357      !!----------------------------------------------------------------------
1358      !!             ***  ROUTINE sbc_cpl_ice_tau  ***
1359      !!
1360      !! ** Purpose :   provide the stress over sea-ice in coupled mode
1361      !!
1362      !! ** Method  :   transform the received stress from the atmosphere into
1363      !!             an atmosphere-ice stress in the (i,j) ocean referencial
1364      !!             and at the velocity point of the sea-ice model (cp_ice_msh):
1365      !!                'C'-grid : i- (j-) components given at U- (V-) point
1366      !!                'I'-grid : B-grid lower-left corner: both components given at I-point
1367      !!
1368      !!                The received stress are :
1369      !!                 - defined by 3 components (if cartesian coordinate)
1370      !!                        or by 2 components (if spherical)
1371      !!                 - oriented along geographical   coordinate (if eastward-northward)
1372      !!                        or  along the local grid coordinate (if local grid)
1373      !!                 - given at U- and V-point, resp.   if received on 2 grids
1374      !!                        or at a same point (T or I) if received on 1 grid
1375      !!                Therefore and if necessary, they are successively
1376      !!             processed in order to obtain them
1377      !!                 first  as  2 components on the sphere
1378      !!                 second as  2 components oriented along the local grid
1379      !!                 third  as  2 components on the cp_ice_msh point
1380      !!
1381      !!                Except in 'oce and ice' case, only one vector stress field
1382      !!             is received. It has already been processed in sbc_cpl_rcv
1383      !!             so that it is now defined as (i,j) components given at U-
1384      !!             and V-points, respectively. Therefore, only the third
1385      !!             transformation is done and only if the ice-grid is a 'I'-grid.
1386      !!
1387      !! ** Action  :   return ptau_i, ptau_j, the stress over the ice at cp_ice_msh point
1388      !!----------------------------------------------------------------------
1389      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_taui   ! i- & j-components of atmos-ice stress [N/m2]
1390      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_tauj   ! at I-point (B-grid) or U & V-point (C-grid)
1391      !!
1392      INTEGER ::   ji, jj   ! dummy loop indices
1393      INTEGER ::   itx      ! index of taux over ice
1394      REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty 
1395      !!----------------------------------------------------------------------
1396      !
1397      IF( nn_timing == 1 )   CALL timing_start('sbc_cpl_ice_tau')
1398      !
1399      CALL wrk_alloc( jpi,jpj,   ztx, zty )
1400
1401      IF( srcv(jpr_itx1)%laction ) THEN   ;   itx =  jpr_itx1   
1402      ELSE                                ;   itx =  jpr_otx1
1403      ENDIF
1404
1405      ! do something only if we just received the stress from atmosphere
1406      IF(  nrcvinfo(itx) == OASIS_Rcv ) THEN
1407         !                                                      ! ======================= !
1408         IF( srcv(jpr_itx1)%laction ) THEN                      !   ice stress received   !
1409            !                                                   ! ======================= !
1410           
1411            IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN            ! 2 components on the sphere
1412               !                                                       ! (cartesian to spherical -> 3 to 2 components)
1413               CALL geo2oce(  frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), frcv(jpr_itz1)%z3(:,:,1),   &
1414                  &          srcv(jpr_itx1)%clgrid, ztx, zty )
1415               frcv(jpr_itx1)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 1st grid
1416               frcv(jpr_ity1)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 1st grid
1417               !
1418               IF( srcv(jpr_itx2)%laction ) THEN
1419                  CALL geo2oce( frcv(jpr_itx2)%z3(:,:,1), frcv(jpr_ity2)%z3(:,:,1), frcv(jpr_itz2)%z3(:,:,1),   &
1420                     &          srcv(jpr_itx2)%clgrid, ztx, zty )
1421                  frcv(jpr_itx2)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 2nd grid
1422                  frcv(jpr_ity2)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 2nd grid
1423               ENDIF
1424               !
1425            ENDIF
1426            !
1427            IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN   ! 2 components oriented along the local grid
1428               !                                                       ! (geographical to local grid -> rotate the components)
1429               CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), srcv(jpr_itx1)%clgrid, 'en->i', ztx )   
1430               IF( srcv(jpr_itx2)%laction ) THEN
1431                  CALL rot_rep( frcv(jpr_itx2)%z3(:,:,1), frcv(jpr_ity2)%z3(:,:,1), srcv(jpr_itx2)%clgrid, 'en->j', zty )   
1432               ELSE
1433                  CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), srcv(jpr_itx1)%clgrid, 'en->j', zty ) 
1434               ENDIF
1435               frcv(jpr_itx1)%z3(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid
1436               frcv(jpr_ity1)%z3(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 1st grid
1437            ENDIF
1438            !                                                   ! ======================= !
1439         ELSE                                                   !     use ocean stress    !
1440            !                                                   ! ======================= !
1441            frcv(jpr_itx1)%z3(:,:,1) = frcv(jpr_otx1)%z3(:,:,1)
1442            frcv(jpr_ity1)%z3(:,:,1) = frcv(jpr_oty1)%z3(:,:,1)
1443            !
1444         ENDIF
1445         !                                                      ! ======================= !
1446         !                                                      !     put on ice grid     !
1447         !                                                      ! ======================= !
1448         !   
1449         !                                                  j+1   j     -----V---F
1450         ! ice stress on ice velocity point (cp_ice_msh)                 !       |
1451         ! (C-grid ==>(U,V) or B-grid ==> I or F)                 j      |   T   U
1452         !                                                               |       |
1453         !                                                   j    j-1   -I-------|
1454         !                                               (for I)         |       |
1455         !                                                              i-1  i   i
1456         !                                                               i      i+1 (for I)
1457         SELECT CASE ( cp_ice_msh )
1458            !
1459         CASE( 'I' )                                         ! B-grid ==> I
1460            SELECT CASE ( srcv(jpr_itx1)%clgrid )
1461            CASE( 'U' )
1462               DO jj = 2, jpjm1                                   ! (U,V) ==> I
1463                  DO ji = 2, jpim1   ! NO vector opt.
1464                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji-1,jj  ,1) + frcv(jpr_itx1)%z3(ji-1,jj-1,1) )
1465                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji  ,jj-1,1) + frcv(jpr_ity1)%z3(ji-1,jj-1,1) )
1466                  END DO
1467               END DO
1468            CASE( 'F' )
1469               DO jj = 2, jpjm1                                   ! F ==> I
1470                  DO ji = 2, jpim1   ! NO vector opt.
1471                     p_taui(ji,jj) = frcv(jpr_itx1)%z3(ji-1,jj-1,1)
1472                     p_tauj(ji,jj) = frcv(jpr_ity1)%z3(ji-1,jj-1,1)
1473                  END DO
1474               END DO
1475            CASE( 'T' )
1476               DO jj = 2, jpjm1                                   ! T ==> I
1477                  DO ji = 2, jpim1   ! NO vector opt.
1478                     p_taui(ji,jj) = 0.25 * ( frcv(jpr_itx1)%z3(ji,jj  ,1) + frcv(jpr_itx1)%z3(ji-1,jj  ,1)   &
1479                        &                   + frcv(jpr_itx1)%z3(ji,jj-1,1) + frcv(jpr_itx1)%z3(ji-1,jj-1,1) ) 
1480                     p_tauj(ji,jj) = 0.25 * ( frcv(jpr_ity1)%z3(ji,jj  ,1) + frcv(jpr_ity1)%z3(ji-1,jj  ,1)   &
1481                        &                   + frcv(jpr_oty1)%z3(ji,jj-1,1) + frcv(jpr_ity1)%z3(ji-1,jj-1,1) )
1482                  END DO
1483               END DO
1484            CASE( 'I' )
1485               p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! I ==> I
1486               p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1487            END SELECT
1488            IF( srcv(jpr_itx1)%clgrid /= 'I' ) THEN
1489               CALL lbc_lnk( p_taui, 'I',  -1. )   ;   CALL lbc_lnk( p_tauj, 'I',  -1. )
1490            ENDIF
1491            !
1492         CASE( 'F' )                                         ! B-grid ==> F
1493            SELECT CASE ( srcv(jpr_itx1)%clgrid )
1494            CASE( 'U' )
1495               DO jj = 2, jpjm1                                   ! (U,V) ==> F
1496                  DO ji = fs_2, fs_jpim1   ! vector opt.
1497                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji,jj,1) + frcv(jpr_itx1)%z3(ji  ,jj+1,1) )
1498                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji,jj,1) + frcv(jpr_ity1)%z3(ji+1,jj  ,1) )
1499                  END DO
1500               END DO
1501            CASE( 'I' )
1502               DO jj = 2, jpjm1                                   ! I ==> F
1503                  DO ji = 2, jpim1   ! NO vector opt.
1504                     p_taui(ji,jj) = frcv(jpr_itx1)%z3(ji+1,jj+1,1)
1505                     p_tauj(ji,jj) = frcv(jpr_ity1)%z3(ji+1,jj+1,1)
1506                  END DO
1507               END DO
1508            CASE( 'T' )
1509               DO jj = 2, jpjm1                                   ! T ==> F
1510                  DO ji = 2, jpim1   ! NO vector opt.
1511                     p_taui(ji,jj) = 0.25 * ( frcv(jpr_itx1)%z3(ji,jj  ,1) + frcv(jpr_itx1)%z3(ji+1,jj  ,1)   &
1512                        &                   + frcv(jpr_itx1)%z3(ji,jj+1,1) + frcv(jpr_itx1)%z3(ji+1,jj+1,1) ) 
1513                     p_tauj(ji,jj) = 0.25 * ( frcv(jpr_ity1)%z3(ji,jj  ,1) + frcv(jpr_ity1)%z3(ji+1,jj  ,1)   &
1514                        &                   + frcv(jpr_ity1)%z3(ji,jj+1,1) + frcv(jpr_ity1)%z3(ji+1,jj+1,1) )
1515                  END DO
1516               END DO
1517            CASE( 'F' )
1518               p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! F ==> F
1519               p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1520            END SELECT
1521            IF( srcv(jpr_itx1)%clgrid /= 'F' ) THEN
1522               CALL lbc_lnk( p_taui, 'F',  -1. )   ;   CALL lbc_lnk( p_tauj, 'F',  -1. )
1523            ENDIF
1524            !
1525         CASE( 'C' )                                         ! C-grid ==> U,V
1526            SELECT CASE ( srcv(jpr_itx1)%clgrid )
1527            CASE( 'U' )
1528               p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! (U,V) ==> (U,V)
1529               p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1530            CASE( 'F' )
1531               DO jj = 2, jpjm1                                   ! F ==> (U,V)
1532                  DO ji = fs_2, fs_jpim1   ! vector opt.
1533                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji,jj,1) + frcv(jpr_itx1)%z3(ji  ,jj-1,1) )
1534                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(jj,jj,1) + frcv(jpr_ity1)%z3(ji-1,jj  ,1) )
1535                  END DO
1536               END DO
1537            CASE( 'T' )
1538               DO jj = 2, jpjm1                                   ! T ==> (U,V)
1539                  DO ji = fs_2, fs_jpim1   ! vector opt.
1540                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj  ,1) + frcv(jpr_itx1)%z3(ji,jj,1) )
1541                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji  ,jj+1,1) + frcv(jpr_ity1)%z3(ji,jj,1) )
1542                  END DO
1543               END DO
1544            CASE( 'I' )
1545               DO jj = 2, jpjm1                                   ! I ==> (U,V)
1546                  DO ji = 2, jpim1   ! NO vector opt.
1547                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj+1,1) + frcv(jpr_itx1)%z3(ji+1,jj  ,1) )
1548                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji+1,jj+1,1) + frcv(jpr_ity1)%z3(ji  ,jj+1,1) )
1549                  END DO
1550               END DO
1551            END SELECT
1552            IF( srcv(jpr_itx1)%clgrid /= 'U' ) THEN
1553               CALL lbc_lnk( p_taui, 'U',  -1. )   ;   CALL lbc_lnk( p_tauj, 'V',  -1. )
1554            ENDIF
1555         END SELECT
1556
1557      ENDIF
1558      !   
1559      CALL wrk_dealloc( jpi,jpj,   ztx, zty )
1560      !
1561      IF( nn_timing == 1 )   CALL timing_stop('sbc_cpl_ice_tau')
1562      !
1563   END SUBROUTINE sbc_cpl_ice_tau
1564   
1565
1566   SUBROUTINE sbc_cpl_ice_flx( p_frld, palbi, psst, pist )
1567      !!----------------------------------------------------------------------
1568      !!             ***  ROUTINE sbc_cpl_ice_flx  ***
1569      !!
1570      !! ** Purpose :   provide the heat and freshwater fluxes of the ocean-ice system
1571      !!
1572      !! ** Method  :   transform the fields received from the atmosphere into
1573      !!             surface heat and fresh water boundary condition for the
1574      !!             ice-ocean system. The following fields are provided:
1575      !!               * total non solar, solar and freshwater fluxes (qns_tot,
1576      !!             qsr_tot and emp_tot) (total means weighted ice-ocean flux)
1577      !!             NB: emp_tot include runoffs and calving.
1578      !!               * fluxes over ice (qns_ice, qsr_ice, emp_ice) where
1579      !!             emp_ice = sublimation - solid precipitation as liquid
1580      !!             precipitation are re-routed directly to the ocean and
1581      !!             calving directly enter the ocean (runoffs are read but included in trasbc.F90)
1582      !!               * solid precipitation (sprecip), used to add to qns_tot
1583      !!             the heat lost associated to melting solid precipitation
1584      !!             over the ocean fraction.
1585      !!               * heat content of rain, snow and evap can also be provided,
1586      !!             otherwise heat flux associated with these mass flux are
1587      !!             guessed (qemp_oce, qemp_ice)
1588      !!
1589      !!             - the fluxes have been separated from the stress as
1590      !!               (a) they are updated at each ice time step compare to
1591      !!               an update at each coupled time step for the stress, and
1592      !!               (b) the conservative computation of the fluxes over the
1593      !!               sea-ice area requires the knowledge of the ice fraction
1594      !!               after the ice advection and before the ice thermodynamics,
1595      !!               so that the stress is updated before the ice dynamics
1596      !!               while the fluxes are updated after it.
1597      !!
1598      !! ** Details
1599      !!             qns_tot = pfrld * qns_oce + ( 1 - pfrld ) * qns_ice   => provided
1600      !!                     + qemp_oce + qemp_ice                         => recalculated and added up to qns
1601      !!
1602      !!             qsr_tot = pfrld * qsr_oce + ( 1 - pfrld ) * qsr_ice   => provided
1603      !!
1604      !!             emp_tot = emp_oce + emp_ice                           => calving is provided and added to emp_tot (and emp_oce).
1605      !!                                                                      runoff (which includes rivers+icebergs) and iceshelf
1606      !!                                                                      are provided but not included in emp here. Only runoff will
1607      !!                                                                      be included in emp in other parts of NEMO code
1608      !! ** Action  :   update at each nf_ice time step:
1609      !!                   qns_tot, qsr_tot  non-solar and solar total heat fluxes
1610      !!                   qns_ice, qsr_ice  non-solar and solar heat fluxes over the ice
1611      !!                   emp_tot           total evaporation - precipitation(liquid and solid) (-calving)
1612      !!                   emp_ice           ice sublimation - solid precipitation over the ice
1613      !!                   dqns_ice          d(non-solar heat flux)/d(Temperature) over the ice
1614      !!                   sprecip           solid precipitation over the ocean 
1615      !!----------------------------------------------------------------------
1616      REAL(wp), INTENT(in   ), DIMENSION(:,:)   ::   p_frld     ! lead fraction                [0 to 1]
1617      ! optional arguments, used only in 'mixed oce-ice' case
1618      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   palbi      ! all skies ice albedo
1619      REAL(wp), INTENT(in   ), DIMENSION(:,:  ), OPTIONAL ::   psst       ! sea surface temperature     [Celsius]
1620      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   pist       ! ice surface temperature     [Kelvin]
1621      !
1622      INTEGER ::   jl         ! dummy loop index
1623      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zcptn, zcptrain, zcptsnw, zicefr, zmsk, zsnw
1624      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice
1625      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice
1626      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice
1627      !!----------------------------------------------------------------------
1628      !
1629      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_flx')
1630      !
1631      CALL wrk_alloc( jpi,jpj,     zcptn, zcptrain, zcptsnw, zicefr, zmsk, zsnw )
1632      CALL wrk_alloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice )
1633      CALL wrk_alloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice )
1634      CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice )
1635
1636      IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
1637      zicefr(:,:) = 1.- p_frld(:,:)
1638      zcptn(:,:) = rcp * sst_m(:,:)
1639      !
1640      !                                                      ! ========================= !
1641      !                                                      !    freshwater budget      !   (emp_tot)
1642      !                                                      ! ========================= !
1643      !
1644      !                                                           ! solid Precipitation                                (sprecip)
1645      !                                                           ! liquid + solid Precipitation                       (tprecip)
1646      !                                                           ! total Evaporation - total Precipitation            (emp_tot)
1647      !                                                           ! sublimation - solid precipitation (cell average)   (emp_ice)
1648      SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
1649      CASE( 'conservative' )   ! received fields: jpr_rain, jpr_snow, jpr_ievp, jpr_tevp
1650         zsprecip(:,:) =   frcv(jpr_snow)%z3(:,:,1)                  ! May need to ensure positive here
1651         ztprecip(:,:) =   frcv(jpr_rain)%z3(:,:,1) + zsprecip(:,:)  ! May need to ensure positive here
1652         zemp_tot(:,:) =   frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:)
1653         zemp_ice(:,:) = ( frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1) ) * zicefr(:,:)
1654      CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
1655         zemp_tot(:,:) = p_frld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + zicefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
1656         zemp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1) * zicefr(:,:)
1657         zsprecip(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_semp)%z3(:,:,1)
1658         ztprecip(:,:) = frcv(jpr_semp)%z3(:,:,1) - frcv(jpr_sbpr)%z3(:,:,1) + zsprecip(:,:)
1659      END SELECT
1660
1661#if defined key_lim3
1662      ! zsnw = snow fraction over ice after wind blowing (=zicefr if no blowing)
1663      zsnw(:,:) = 0._wp  ;  CALL lim_thd_snwblow( p_frld, zsnw )
1664     
1665      ! --- evaporation minus precipitation corrected (because of wind blowing on snow) --- !
1666      zemp_ice(:,:) = zemp_ice(:,:) + zsprecip(:,:) * ( zicefr(:,:) - zsnw(:,:) )  ! emp_ice = A * sublimation - zsnw * sprecip
1667      zemp_oce(:,:) = zemp_tot(:,:) - zemp_ice(:,:)                                ! emp_oce = emp_tot - emp_ice
1668
1669      ! --- evaporation over ocean (used later for qemp) --- !
1670      zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:)
1671
1672      ! --- evaporation over ice (kg/m2/s) --- !
1673      zevap_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1)
1674      ! since the sensitivity of evap to temperature (devap/dT) is not prescribed by the atmosphere, we set it to 0
1675      ! therefore, sublimation is not redistributed over the ice categories when no subgrid scale fluxes are provided by atm.
1676      zdevap_ice(:,:) = 0._wp
1677     
1678      ! --- Continental fluxes --- !
1679      IF( srcv(jpr_rnf)%laction ) THEN   ! runoffs (included in emp later on)
1680         rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
1681      ENDIF
1682      IF( srcv(jpr_cal)%laction ) THEN   ! calving (put in emp_tot and emp_oce)
1683         zemp_tot(:,:) = zemp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
1684         zemp_oce(:,:) = zemp_oce(:,:) - frcv(jpr_cal)%z3(:,:,1)
1685      ENDIF
1686      IF( srcv(jpr_icb)%laction ) THEN   ! iceberg added to runoffs
1687         fwficb(:,:) = frcv(jpr_icb)%z3(:,:,1)
1688         rnf(:,:)    = rnf(:,:) + fwficb(:,:)
1689      ENDIF
1690      IF( srcv(jpr_isf)%laction ) THEN   ! iceshelf (fwfisf <0 mean melting)
1691        fwfisf(:,:) = - frcv(jpr_isf)%z3(:,:,1) 
1692      ENDIF
1693
1694      IF( ln_mixcpl ) THEN
1695         emp_tot(:,:) = emp_tot(:,:) * xcplmask(:,:,0) + zemp_tot(:,:) * zmsk(:,:)
1696         emp_ice(:,:) = emp_ice(:,:) * xcplmask(:,:,0) + zemp_ice(:,:) * zmsk(:,:)
1697         emp_oce(:,:) = emp_oce(:,:) * xcplmask(:,:,0) + zemp_oce(:,:) * zmsk(:,:)
1698         sprecip(:,:) = sprecip(:,:) * xcplmask(:,:,0) + zsprecip(:,:) * zmsk(:,:)
1699         tprecip(:,:) = tprecip(:,:) * xcplmask(:,:,0) + ztprecip(:,:) * zmsk(:,:)
1700         DO jl=1,jpl
1701            evap_ice (:,:,jl) = evap_ice (:,:,jl) * xcplmask(:,:,0) + zevap_ice (:,:) * zmsk(:,:)
1702            devap_ice(:,:,jl) = devap_ice(:,:,jl) * xcplmask(:,:,0) + zdevap_ice(:,:) * zmsk(:,:)
1703         ENDDO
1704      ELSE
1705         emp_tot(:,:) =         zemp_tot(:,:)
1706         emp_ice(:,:) =         zemp_ice(:,:)
1707         emp_oce(:,:) =         zemp_oce(:,:)     
1708         sprecip(:,:) =         zsprecip(:,:)
1709         tprecip(:,:) =         ztprecip(:,:)
1710         DO jl=1,jpl
1711            evap_ice (:,:,jl) = zevap_ice (:,:)
1712            devap_ice(:,:,jl) = zdevap_ice(:,:)
1713         ENDDO
1714      ENDIF
1715
1716#else
1717      zsnw(:,:) = zicefr(:,:)
1718      ! --- Continental fluxes --- !
1719      IF( srcv(jpr_rnf)%laction ) THEN   ! runoffs (included in emp later on)
1720         rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
1721      ENDIF
1722      IF( srcv(jpr_cal)%laction ) THEN   ! calving (put in emp_tot)
1723         zemp_tot(:,:) = zemp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
1724      ENDIF
1725      IF( srcv(jpr_icb)%laction ) THEN   ! iceberg added to runoffs
1726         fwficb(:,:) = frcv(jpr_icb)%z3(:,:,1)
1727         rnf(:,:)    = rnf(:,:) + fwficb(:,:)
1728      ENDIF
1729      IF( srcv(jpr_isf)%laction ) THEN   ! iceshelf (fwfisf <0 mean melting)
1730        fwfisf(:,:) = - frcv(jpr_isf)%z3(:,:,1)
1731      ENDIF
1732
1733      IF( ln_mixcpl ) THEN
1734         emp_tot(:,:) = emp_tot(:,:) * xcplmask(:,:,0) + zemp_tot(:,:) * zmsk(:,:)
1735         emp_ice(:,:) = emp_ice(:,:) * xcplmask(:,:,0) + zemp_ice(:,:) * zmsk(:,:)
1736         sprecip(:,:) = sprecip(:,:) * xcplmask(:,:,0) + zsprecip(:,:) * zmsk(:,:)
1737         tprecip(:,:) = tprecip(:,:) * xcplmask(:,:,0) + ztprecip(:,:) * zmsk(:,:)
1738      ELSE
1739         emp_tot(:,:) =                                  zemp_tot(:,:)
1740         emp_ice(:,:) =                                  zemp_ice(:,:)
1741         sprecip(:,:) =                                  zsprecip(:,:)
1742         tprecip(:,:) =                                  ztprecip(:,:)
1743      ENDIF
1744
1745#endif
1746      ! outputs
1747!!      IF( srcv(jpr_rnf)%laction )   CALL iom_put( 'runoffs' , rnf(:,:) * tmask(:,:,1)                                 )  ! runoff
1748!!      IF( srcv(jpr_isf)%laction )   CALL iom_put( 'iceshelf_cea', -fwfisf(:,:) * tmask(:,:,1)                         )  ! iceshelf
1749      IF( srcv(jpr_cal)%laction )   CALL iom_put( 'calving_cea' , frcv(jpr_cal)%z3(:,:,1) * tmask(:,:,1)                )  ! calving
1750      IF( srcv(jpr_icb)%laction )   CALL iom_put( 'iceberg_cea' , frcv(jpr_icb)%z3(:,:,1) * tmask(:,:,1)                )  ! icebergs
1751      IF( iom_use('snowpre') )      CALL iom_put( 'snowpre'     , sprecip(:,:)                                          )  ! Snow
1752      IF( iom_use('precip') )       CALL iom_put( 'precip'      , tprecip(:,:)                                          )  ! total  precipitation
1753      IF( iom_use('rain') )         CALL iom_put( 'rain'        , tprecip(:,:) - sprecip(:,:)                           )  ! liquid precipitation
1754      IF( iom_use('snow_ao_cea') )  CALL iom_put( 'snow_ao_cea' , sprecip(:,:) * ( 1._wp - zsnw(:,:) )                  )  ! Snow over ice-free ocean  (cell average)
1755      IF( iom_use('snow_ai_cea') )  CALL iom_put( 'snow_ai_cea' , sprecip(:,:) *           zsnw(:,:)                    )  ! Snow over sea-ice         (cell average)
1756      IF( iom_use('subl_ai_cea') )  CALL iom_put( 'subl_ai_cea' , frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) * tmask(:,:,1) )  ! Sublimation over sea-ice (cell average)
1757      IF( iom_use('evap_ao_cea') )  CALL iom_put( 'evap_ao_cea' , ( frcv(jpr_tevp)%z3(:,:,1)  &
1758         &                                                        - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) ) * tmask(:,:,1) )  ! ice-free oce evap (cell average)
1759      ! note: runoff output is done in sbcrnf (which includes icebergs too) and iceshelf output is done in sbcisf
1760      !
1761      !                                                      ! ========================= !
1762      SELECT CASE( TRIM( sn_rcv_qns%cldes ) )                !   non solar heat fluxes   !   (qns)
1763      !                                                      ! ========================= !
1764      CASE( 'oce only' )         ! the required field is directly provided
1765         zqns_tot(:,:) = frcv(jpr_qnsoce)%z3(:,:,1)
1766      CASE( 'conservative' )     ! the required fields are directly provided
1767         zqns_tot(:,:) = frcv(jpr_qnsmix)%z3(:,:,1)
1768         IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
1769            zqns_ice(:,:,1:jpl) = frcv(jpr_qnsice)%z3(:,:,1:jpl)
1770         ELSE
1771            DO jl=1,jpl
1772               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1) ! Set all category values equal
1773            ENDDO
1774         ENDIF
1775      CASE( 'oce and ice' )      ! the total flux is computed from ocean and ice fluxes
1776         zqns_tot(:,:) =  p_frld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1)
1777         IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
1778            DO jl=1,jpl
1779               zqns_tot(:,:   ) = zqns_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qnsice)%z3(:,:,jl)   
1780               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,jl)
1781            ENDDO
1782         ELSE
1783            qns_tot(:,:) = qns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
1784            DO jl=1,jpl
1785               zqns_tot(:,:   ) = zqns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
1786               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
1787            ENDDO
1788         ENDIF
1789      CASE( 'mixed oce-ice' )    ! the ice flux is cumputed from the total flux, the SST and ice informations
1790! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
1791         zqns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
1792         zqns_ice(:,:,1) = frcv(jpr_qnsmix)%z3(:,:,1)    &
1793            &            + frcv(jpr_dqnsdt)%z3(:,:,1) * ( pist(:,:,1) - ( (rt0 + psst(:,:  ) ) * p_frld(:,:)   &
1794            &                                           + pist(:,:,1) * zicefr(:,:) ) )
1795      END SELECT
1796      !                                     
1797      ! --- calving (removed from qns_tot) --- !
1798      IF( srcv(jpr_cal)%laction )   zqns_tot(:,:) = zqns_tot(:,:) - frcv(jpr_cal)%z3(:,:,1) * lfus  ! remove latent heat of calving
1799                                                                                                    ! we suppose it melts at 0deg, though it should be temp. of surrounding ocean
1800      ! --- iceberg (removed from qns_tot) --- !
1801      IF( srcv(jpr_icb)%laction )   zqns_tot(:,:) = zqns_tot(:,:) - frcv(jpr_icb)%z3(:,:,1) * lfus  ! remove latent heat of iceberg melting
1802
1803#if defined key_lim3     
1804      ! --- non solar flux over ocean --- !
1805      !         note: p_frld cannot be = 0 since we limit the ice concentration to amax
1806      zqns_oce = 0._wp
1807      WHERE( p_frld /= 0._wp )  zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / p_frld(:,:)
1808
1809      ! Heat content per unit mass of snow (J/kg)
1810      WHERE( SUM( a_i, dim=3 ) > 1.e-10 )   ;   zcptsnw(:,:) = cpic * SUM( (tn_ice - rt0) * a_i, dim=3 ) / SUM( a_i, dim=3 )
1811      ELSEWHERE                             ;   zcptsnw(:,:) = zcptn(:,:)
1812      ENDWHERE
1813      ! Heat content per unit mass of rain (J/kg)
1814      zcptrain(:,:) = rcp * ( SUM( (tn_ice(:,:,:) - rt0) * a_i(:,:,:), dim=3 ) + sst_m(:,:) * p_frld(:,:) ) 
1815
1816      ! --- enthalpy of snow precip over ice in J/m3 (to be used in 1D-thermo) --- !
1817      zqprec_ice(:,:) = rhosn * ( zcptsnw(:,:) - lfus )
1818
1819      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- !
1820      DO jl = 1, jpl
1821         zqevap_ice(:,:,jl) = 0._wp ! should be -evap * ( ( Tice - rt0 ) * cpic ) but atm. does not take it into account
1822      END DO
1823
1824      ! --- heat flux associated with emp (W/m2) --- !
1825      zqemp_oce(:,:) = -  zevap_oce(:,:)                                      *   zcptn   (:,:)   &        ! evap
1826         &             + ( ztprecip(:,:) - zsprecip(:,:) )                    *   zcptrain(:,:)   &        ! liquid precip
1827         &             +   zsprecip(:,:)                   * ( 1._wp - zsnw ) * ( zcptsnw (:,:) - lfus )   ! solid precip over ocean + snow melting
1828      zqemp_ice(:,:) =     zsprecip(:,:)                   * zsnw             * ( zcptsnw (:,:) - lfus )   ! solid precip over ice (qevap_ice=0 since atm. does not take it into account)
1829!!    zqemp_ice(:,:) = -   frcv(jpr_ievp)%z3(:,:,1)        * zicefr(:,:)      *   zcptsnw (:,:)   &        ! ice evap
1830!!       &             +   zsprecip(:,:)                   * zsnw             * zqprec_ice(:,:) * r1_rhosn ! solid precip over ice
1831     
1832      ! --- total non solar flux (including evap/precip) --- !
1833      zqns_tot(:,:) = zqns_tot(:,:) + zqemp_ice(:,:) + zqemp_oce(:,:)
1834
1835      ! --- in case both coupled/forced are active, we must mix values --- !
1836      IF( ln_mixcpl ) THEN
1837         qns_tot(:,:) = qns_tot(:,:) * xcplmask(:,:,0) + zqns_tot(:,:)* zmsk(:,:)
1838         qns_oce(:,:) = qns_oce(:,:) * xcplmask(:,:,0) + zqns_oce(:,:)* zmsk(:,:)
1839         DO jl=1,jpl
1840            qns_ice  (:,:,jl) = qns_ice  (:,:,jl) * xcplmask(:,:,0) +  zqns_ice  (:,:,jl)* zmsk(:,:)
1841            qevap_ice(:,:,jl) = qevap_ice(:,:,jl) * xcplmask(:,:,0) +  zqevap_ice(:,:,jl)* zmsk(:,:)
1842         ENDDO
1843         qprec_ice(:,:) = qprec_ice(:,:) * xcplmask(:,:,0) + zqprec_ice(:,:)* zmsk(:,:)
1844         qemp_oce (:,:) =  qemp_oce(:,:) * xcplmask(:,:,0) +  zqemp_oce(:,:)* zmsk(:,:)
1845         qemp_ice (:,:) =  qemp_ice(:,:) * xcplmask(:,:,0) +  zqemp_ice(:,:)* zmsk(:,:)
1846      ELSE
1847         qns_tot  (:,:  ) = zqns_tot  (:,:  )
1848         qns_oce  (:,:  ) = zqns_oce  (:,:  )
1849         qns_ice  (:,:,:) = zqns_ice  (:,:,:)
1850         qevap_ice(:,:,:) = zqevap_ice(:,:,:)
1851         qprec_ice(:,:  ) = zqprec_ice(:,:  )
1852         qemp_oce (:,:  ) = zqemp_oce (:,:  )
1853         qemp_ice (:,:  ) = zqemp_ice (:,:  )
1854      ENDIF
1855
1856#else
1857      zcptsnw (:,:) = zcptn(:,:)
1858      zcptrain(:,:) = zcptn(:,:)
1859     
1860      ! clem: this formulation is certainly wrong... but better than it was...
1861      zqns_tot(:,:) = zqns_tot(:,:)                            &          ! zqns_tot update over free ocean with:
1862         &          - (  p_frld(:,:) * zsprecip(:,:) * lfus )  &          ! remove the latent heat flux of solid precip. melting
1863         &          - (  zemp_tot(:,:)                         &          ! remove the heat content of mass flux (assumed to be at SST)
1864         &             - zemp_ice(:,:) ) * zcptn(:,:) 
1865
1866     IF( ln_mixcpl ) THEN
1867         qns_tot(:,:) = qns(:,:) * p_frld(:,:) + SUM( qns_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk
1868         qns_tot(:,:) = qns_tot(:,:) * xcplmask(:,:,0) +  zqns_tot(:,:)* zmsk(:,:)
1869         DO jl=1,jpl
1870            qns_ice(:,:,jl) = qns_ice(:,:,jl) * xcplmask(:,:,0) +  zqns_ice(:,:,jl)* zmsk(:,:)
1871         ENDDO
1872      ELSE
1873         qns_tot(:,:  ) = zqns_tot(:,:  )
1874         qns_ice(:,:,:) = zqns_ice(:,:,:)
1875      ENDIF
1876
1877#endif
1878      ! outputs
1879      IF( srcv(jpr_cal)%laction )    CALL iom_put('hflx_cal_cea' ,   &
1880                                             &   - frcv(jpr_cal)%z3(:,:,1) * lfus) ! latent heat from calving
1881      IF( srcv(jpr_icb)%laction )    CALL iom_put('hflx_icb_cea' , - frcv(jpr_icb)%z3(:,:,1) * lfus) ! latent heat from icebergs melting
1882      IF( iom_use('hflx_snow_cea') ) CALL iom_put('hflx_snow_cea',  sprecip(:,:) * ( zcptsnw(:,:) - Lfus )) ! heat flux from snow (cell average)
1883      IF( iom_use('hflx_rain_cea') ) CALL iom_put('hflx_rain_cea',( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:)) ! heat flux from rain (cell average)
1884      IF( iom_use('hflx_evap_cea') ) CALL iom_put('hflx_evap_cea',(frcv(jpr_tevp)%z3(:,:,1)   &
1885                                             &    - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) & ! heat flux from from evap (cell average)
1886         &                                                        ) * zcptn(:,:) * tmask(:,:,1) )
1887      IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea',sprecip(:,:)   &
1888                                             & * (zcptsnw(:,:) - Lfus) * (1._wp - zsnw(:,:))   ) ! heat flux from snow (over ocean)
1889      IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea',sprecip(:,:) * (zcptsnw(:,:) - Lfus)   &
1890                                             & *          zsnw(:,:)    ) ! heat flux from snow (over ice)
1891      ! note: hflx for runoff and iceshelf are done in sbcrnf and sbcisf resp.
1892      !
1893      !                                                      ! ========================= !
1894      SELECT CASE( TRIM( sn_rcv_qsr%cldes ) )                !      solar heat fluxes    !   (qsr)
1895      !                                                      ! ========================= !
1896      CASE( 'oce only' )
1897         zqsr_tot(:,:  ) = MAX( 0._wp , frcv(jpr_qsroce)%z3(:,:,1) )
1898      CASE( 'conservative' )
1899         zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
1900         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
1901            zqsr_ice(:,:,1:jpl) = frcv(jpr_qsrice)%z3(:,:,1:jpl)
1902         ELSE
1903            ! Set all category values equal for the moment
1904            DO jl=1,jpl
1905               zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
1906            ENDDO
1907         ENDIF
1908         zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
1909         zqsr_ice(:,:,1) = frcv(jpr_qsrice)%z3(:,:,1)
1910      CASE( 'oce and ice' )
1911         zqsr_tot(:,:  ) =  p_frld(:,:) * frcv(jpr_qsroce)%z3(:,:,1)
1912         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
1913            DO jl=1,jpl
1914               zqsr_tot(:,:   ) = zqsr_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qsrice)%z3(:,:,jl)   
1915               zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,jl)
1916            ENDDO
1917         ELSE
1918            qsr_tot(:,:   ) = qsr_tot(:,:) + zicefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
1919            DO jl=1,jpl
1920               zqsr_tot(:,:   ) = zqsr_tot(:,:) + zicefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
1921               zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
1922            ENDDO
1923         ENDIF
1924      CASE( 'mixed oce-ice' )
1925         zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
1926! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
1927!       Create solar heat flux over ice using incoming solar heat flux and albedos
1928!       ( see OASIS3 user guide, 5th edition, p39 )
1929         zqsr_ice(:,:,1) = frcv(jpr_qsrmix)%z3(:,:,1) * ( 1.- palbi(:,:,1) )   &
1930            &            / (  1.- ( albedo_oce_mix(:,:  ) * p_frld(:,:)       &
1931            &                     + palbi         (:,:,1) * zicefr(:,:) ) )
1932      END SELECT
1933      IF( ln_dm2dc .AND. ln_cpl ) THEN   ! modify qsr to include the diurnal cycle
1934         zqsr_tot(:,:  ) = sbc_dcy( zqsr_tot(:,:  ) )
1935         DO jl=1,jpl
1936            zqsr_ice(:,:,jl) = sbc_dcy( zqsr_ice(:,:,jl) )
1937         ENDDO
1938      ENDIF
1939
1940#if defined key_lim3
1941      ! --- solar flux over ocean --- !
1942      !         note: p_frld cannot be = 0 since we limit the ice concentration to amax
1943      zqsr_oce = 0._wp
1944      WHERE( p_frld /= 0._wp )  zqsr_oce(:,:) = ( zqsr_tot(:,:) - SUM( a_i * zqsr_ice, dim=3 ) ) / p_frld(:,:)
1945
1946      IF( ln_mixcpl ) THEN   ;   qsr_oce(:,:) = qsr_oce(:,:) * xcplmask(:,:,0) +  zqsr_oce(:,:)* zmsk(:,:)
1947      ELSE                   ;   qsr_oce(:,:) = zqsr_oce(:,:)   ;   ENDIF
1948#endif
1949
1950      IF( ln_mixcpl ) THEN
1951         qsr_tot(:,:) = qsr(:,:) * p_frld(:,:) + SUM( qsr_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk
1952         qsr_tot(:,:) = qsr_tot(:,:) * xcplmask(:,:,0) +  zqsr_tot(:,:)* zmsk(:,:)
1953         DO jl=1,jpl
1954            qsr_ice(:,:,jl) = qsr_ice(:,:,jl) * xcplmask(:,:,0) +  zqsr_ice(:,:,jl)* zmsk(:,:)
1955         ENDDO
1956      ELSE
1957         qsr_tot(:,:  ) = zqsr_tot(:,:  )
1958         qsr_ice(:,:,:) = zqsr_ice(:,:,:)
1959      ENDIF
1960
1961      !                                                      ! ========================= !
1962      SELECT CASE( TRIM( sn_rcv_dqnsdt%cldes ) )             !          d(qns)/dt        !
1963      !                                                      ! ========================= !
1964      CASE ('coupled')
1965         IF ( TRIM(sn_rcv_dqnsdt%clcat) == 'yes' ) THEN
1966            zdqns_ice(:,:,1:jpl) = frcv(jpr_dqnsdt)%z3(:,:,1:jpl)
1967         ELSE
1968            ! Set all category values equal for the moment
1969            DO jl=1,jpl
1970               zdqns_ice(:,:,jl) = frcv(jpr_dqnsdt)%z3(:,:,1)
1971            ENDDO
1972         ENDIF
1973      END SELECT
1974     
1975      IF( ln_mixcpl ) THEN
1976         DO jl=1,jpl
1977            dqns_ice(:,:,jl) = dqns_ice(:,:,jl) * xcplmask(:,:,0) + zdqns_ice(:,:,jl) * zmsk(:,:)
1978         ENDDO
1979      ELSE
1980         dqns_ice(:,:,:) = zdqns_ice(:,:,:)
1981      ENDIF
1982     
1983      !                                                      ! ========================= !
1984      SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) )             !    topmelt and botmelt    !
1985      !                                                      ! ========================= !
1986      CASE ('coupled')
1987         topmelt(:,:,:)=frcv(jpr_topm)%z3(:,:,:)
1988         botmelt(:,:,:)=frcv(jpr_botm)%z3(:,:,:)
1989      END SELECT
1990
1991      ! Surface transimission parameter io (Maykut Untersteiner , 1971 ; Ebert and Curry, 1993 )
1992      ! Used for LIM2 and LIM3
1993      ! Coupled case: since cloud cover is not received from atmosphere
1994      !               ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
1995      fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
1996      fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
1997
1998      CALL wrk_dealloc( jpi,jpj,     zcptn, zcptrain, zcptsnw, zicefr, zmsk, zsnw )
1999      CALL wrk_dealloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice )
2000      CALL wrk_dealloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice )
2001      CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice )
2002      !
2003      IF( nn_timing == 1 )   CALL timing_stop('sbc_cpl_ice_flx')
2004      !
2005   END SUBROUTINE sbc_cpl_ice_flx
2006   
2007   
2008   SUBROUTINE sbc_cpl_snd( kt )
2009      !!----------------------------------------------------------------------
2010      !!             ***  ROUTINE sbc_cpl_snd  ***
2011      !!
2012      !! ** Purpose :   provide the ocean-ice informations to the atmosphere
2013      !!
2014      !! ** Method  :   send to the atmosphere through a call to cpl_snd
2015      !!              all the needed fields (as defined in sbc_cpl_init)
2016      !!----------------------------------------------------------------------
2017      INTEGER, INTENT(in) ::   kt
2018      !
2019      INTEGER ::   ji, jj, jl   ! dummy loop indices
2020      INTEGER ::   isec, info   ! local integer
2021      REAL(wp) ::   zumax, zvmax
2022      REAL(wp), POINTER, DIMENSION(:,:)   ::   zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1
2023      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmp3, ztmp4   
2024      !!----------------------------------------------------------------------
2025      !
2026      IF( nn_timing == 1 )   CALL timing_start('sbc_cpl_snd')
2027      !
2028      CALL wrk_alloc( jpi,jpj,       zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
2029      CALL wrk_alloc( jpi,jpj,jpl,   ztmp3, ztmp4 )
2030
2031      isec = ( kt - nit000 ) * NINT( rdt )        ! date of exchanges
2032
2033      zfr_l(:,:) = 1.- fr_i(:,:)
2034      !                                                      ! ------------------------- !
2035      !                                                      !    Surface temperature    !   in Kelvin
2036      !                                                      ! ------------------------- !
2037      IF( ssnd(jps_toce)%laction .OR. ssnd(jps_tice)%laction .OR. ssnd(jps_tmix)%laction ) THEN
2038         
2039         IF ( nn_components == jp_iam_opa ) THEN
2040            ztmp1(:,:) = tsn(:,:,1,jp_tem)   ! send temperature as it is (potential or conservative) -> use of l_useCT on the received part
2041         ELSE
2042            ! we must send the surface potential temperature
2043            IF( l_useCT )  THEN    ;   ztmp1(:,:) = eos_pt_from_ct( tsn(:,:,1,jp_tem), tsn(:,:,1,jp_sal) )
2044            ELSE                    ;   ztmp1(:,:) = tsn(:,:,1,jp_tem)
2045            ENDIF
2046            !
2047            SELECT CASE( sn_snd_temp%cldes)
2048            CASE( 'oce only'             )   ;   ztmp1(:,:) =   ztmp1(:,:) + rt0
2049            CASE( 'oce and ice'          )   ;   ztmp1(:,:) =   ztmp1(:,:) + rt0
2050               SELECT CASE( sn_snd_temp%clcat )
2051               CASE( 'yes' )   
2052                  ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl)
2053               CASE( 'no' )
2054                  WHERE( SUM( a_i, dim=3 ) /= 0. )
2055                     ztmp3(:,:,1) = SUM( tn_ice * a_i, dim=3 ) / SUM( a_i, dim=3 )
2056                  ELSEWHERE
2057                     ztmp3(:,:,1) = rt0
2058                  END WHERE
2059               CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
2060               END SELECT
2061            CASE( 'weighted oce and ice' )   ;   ztmp1(:,:) = ( ztmp1(:,:) + rt0 ) * zfr_l(:,:)   
2062               SELECT CASE( sn_snd_temp%clcat )
2063               CASE( 'yes' )   
2064                  ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
2065               CASE( 'no' )
2066                  ztmp3(:,:,:) = 0.0
2067                  DO jl=1,jpl
2068                     ztmp3(:,:,1) = ztmp3(:,:,1) + tn_ice(:,:,jl) * a_i(:,:,jl)
2069                  ENDDO
2070               CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
2071               END SELECT
2072            CASE( 'mixed oce-ice'        )   
2073               ztmp1(:,:) = ( ztmp1(:,:) + rt0 ) * zfr_l(:,:) 
2074               DO jl=1,jpl
2075                  ztmp1(:,:) = ztmp1(:,:) + tn_ice(:,:,jl) * a_i(:,:,jl)
2076               ENDDO
2077            CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%cldes' )
2078            END SELECT
2079         ENDIF
2080         IF( ssnd(jps_toce)%laction )   CALL cpl_snd( jps_toce, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
2081         IF( ssnd(jps_tice)%laction )   CALL cpl_snd( jps_tice, isec, ztmp3, info )
2082         IF( ssnd(jps_tmix)%laction )   CALL cpl_snd( jps_tmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
2083      ENDIF
2084      !                                                      ! ------------------------- !
2085      !                                                      !           Albedo          !
2086      !                                                      ! ------------------------- !
2087      IF( ssnd(jps_albice)%laction ) THEN                         ! ice
2088          SELECT CASE( sn_snd_alb%cldes )
2089          CASE( 'ice' )
2090             SELECT CASE( sn_snd_alb%clcat )
2091             CASE( 'yes' )   
2092                ztmp3(:,:,1:jpl) = alb_ice(:,:,1:jpl)
2093             CASE( 'no' )
2094                WHERE( SUM( a_i, dim=3 ) /= 0. )
2095                   ztmp1(:,:) = SUM( alb_ice (:,:,1:jpl) * a_i(:,:,1:jpl), dim=3 ) / SUM( a_i(:,:,1:jpl), dim=3 )
2096                ELSEWHERE
2097                   ztmp1(:,:) = albedo_oce_mix(:,:)
2098                END WHERE
2099             CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_alb%clcat' )
2100             END SELECT
2101          CASE( 'weighted ice' )   ;
2102             SELECT CASE( sn_snd_alb%clcat )
2103             CASE( 'yes' )   
2104                ztmp3(:,:,1:jpl) =  alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
2105             CASE( 'no' )
2106                WHERE( fr_i (:,:) > 0. )
2107                   ztmp1(:,:) = SUM (  alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl), dim=3 )
2108                ELSEWHERE
2109                   ztmp1(:,:) = 0.
2110                END WHERE
2111             CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_ice%clcat' )
2112             END SELECT
2113          CASE default      ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_alb%cldes' )
2114         END SELECT
2115
2116         SELECT CASE( sn_snd_alb%clcat )
2117            CASE( 'yes' )   
2118               CALL cpl_snd( jps_albice, isec, ztmp3, info )      !-> MV this has never been checked in coupled mode
2119            CASE( 'no'  )   
2120               CALL cpl_snd( jps_albice, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info ) 
2121         END SELECT
2122      ENDIF
2123
2124      IF( ssnd(jps_albmix)%laction ) THEN                         ! mixed ice-ocean
2125         ztmp1(:,:) = albedo_oce_mix(:,:) * zfr_l(:,:)
2126         DO jl=1,jpl
2127            ztmp1(:,:) = ztmp1(:,:) + alb_ice(:,:,jl) * a_i(:,:,jl)
2128         ENDDO
2129         CALL cpl_snd( jps_albmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
2130      ENDIF
2131      !                                                      ! ------------------------- !
2132      !                                                      !  Ice fraction & Thickness !
2133      !                                                      ! ------------------------- !
2134      ! Send ice fraction field to atmosphere
2135      IF( ssnd(jps_fice)%laction ) THEN
2136         SELECT CASE( sn_snd_thick%clcat )
2137         CASE( 'yes' )   ;   ztmp3(:,:,1:jpl) =  a_i(:,:,1:jpl)
2138         CASE( 'no'  )   ;   ztmp3(:,:,1    ) = fr_i(:,:      )
2139         CASE default    ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
2140         END SELECT
2141         IF( ssnd(jps_fice)%laction )   CALL cpl_snd( jps_fice, isec, ztmp3, info )
2142      ENDIF
2143     
2144      ! Send ice fraction field to OPA (sent by SAS in SAS-OPA coupling)
2145      IF( ssnd(jps_fice2)%laction ) THEN
2146         ztmp3(:,:,1) = fr_i(:,:)
2147         IF( ssnd(jps_fice2)%laction )   CALL cpl_snd( jps_fice2, isec, ztmp3, info )
2148      ENDIF
2149
2150      ! Send ice and snow thickness field
2151      IF( ssnd(jps_hice)%laction .OR. ssnd(jps_hsnw)%laction ) THEN
2152         SELECT CASE( sn_snd_thick%cldes)
2153         CASE( 'none'                  )       ! nothing to do
2154         CASE( 'weighted ice and snow' )   
2155            SELECT CASE( sn_snd_thick%clcat )
2156            CASE( 'yes' )   
2157               ztmp3(:,:,1:jpl) =  ht_i(:,:,1:jpl) * a_i(:,:,1:jpl)
2158               ztmp4(:,:,1:jpl) =  ht_s(:,:,1:jpl) * a_i(:,:,1:jpl)
2159            CASE( 'no' )
2160               ztmp3(:,:,:) = 0.0   ;  ztmp4(:,:,:) = 0.0
2161               DO jl=1,jpl
2162                  ztmp3(:,:,1) = ztmp3(:,:,1) + ht_i(:,:,jl) * a_i(:,:,jl)
2163                  ztmp4(:,:,1) = ztmp4(:,:,1) + ht_s(:,:,jl) * a_i(:,:,jl)
2164               ENDDO
2165            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
2166            END SELECT
2167         CASE( 'ice and snow'         )   
2168            SELECT CASE( sn_snd_thick%clcat )
2169            CASE( 'yes' )
2170               ztmp3(:,:,1:jpl) = ht_i(:,:,1:jpl)
2171               ztmp4(:,:,1:jpl) = ht_s(:,:,1:jpl)
2172            CASE( 'no' )
2173               WHERE( SUM( a_i, dim=3 ) /= 0. )
2174                  ztmp3(:,:,1) = SUM( ht_i * a_i, dim=3 ) / SUM( a_i, dim=3 )
2175                  ztmp4(:,:,1) = SUM( ht_s * a_i, dim=3 ) / SUM( a_i, dim=3 )
2176               ELSEWHERE
2177                 ztmp3(:,:,1) = 0.
2178                 ztmp4(:,:,1) = 0.
2179               END WHERE
2180            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
2181            END SELECT
2182         CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%cldes' )
2183         END SELECT
2184         IF( ssnd(jps_hice)%laction )   CALL cpl_snd( jps_hice, isec, ztmp3, info )
2185         IF( ssnd(jps_hsnw)%laction )   CALL cpl_snd( jps_hsnw, isec, ztmp4, info )
2186      ENDIF
2187      !                                                      ! ------------------------- !
2188      !                                                      !  CO2 flux from PISCES     !
2189      !                                                      ! ------------------------- !
2190      IF( ssnd(jps_co2)%laction .AND. l_co2cpl )   CALL cpl_snd( jps_co2, isec, RESHAPE ( oce_co2, (/jpi,jpj,1/) ) , info )
2191      !
2192      !                                                      ! ------------------------- !
2193      IF( ssnd(jps_ocx1)%laction ) THEN                      !      Surface current      !
2194         !                                                   ! ------------------------- !
2195         !   
2196         !                                                  j+1   j     -----V---F
2197         ! surface velocity always sent from T point                     !       |
2198         !                                                        j      |   T   U
2199         !                                                               |       |
2200         !                                                   j    j-1   -I-------|
2201         !                                               (for I)         |       |
2202         !                                                              i-1  i   i
2203         !                                                               i      i+1 (for I)
2204         IF( nn_components == jp_iam_opa ) THEN
2205            zotx1(:,:) = un(:,:,1) 
2206            zoty1(:,:) = vn(:,:,1) 
2207         ELSE       
2208            SELECT CASE( TRIM( sn_snd_crt%cldes ) )
2209            CASE( 'oce only'             )      ! C-grid ==> T
2210               DO jj = 2, jpjm1
2211                  DO ji = fs_2, fs_jpim1   ! vector opt.
2212                     zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj  ,1) )
2213                     zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji  ,jj-1,1) ) 
2214                  END DO
2215               END DO
2216            CASE( 'weighted oce and ice' )   
2217               SELECT CASE ( cp_ice_msh )
2218               CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
2219                  DO jj = 2, jpjm1
2220                     DO ji = fs_2, fs_jpim1   ! vector opt.
2221                        zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj) 
2222                        zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)
2223                        zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
2224                        zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
2225                     END DO
2226                  END DO
2227               CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
2228                  DO jj = 2, jpjm1
2229                     DO ji = 2, jpim1   ! NO vector opt.
2230                        zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj) 
2231                        zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj) 
2232                        zitx1(ji,jj) = 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &
2233                           &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
2234                        zity1(ji,jj) = 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &
2235                           &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
2236                     END DO
2237                  END DO
2238               CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
2239                  DO jj = 2, jpjm1
2240                     DO ji = 2, jpim1   ! NO vector opt.
2241                        zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj) 
2242                        zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj) 
2243                        zitx1(ji,jj) = 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &
2244                           &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
2245                        zity1(ji,jj) = 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &
2246                           &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
2247                     END DO
2248                  END DO
2249               END SELECT
2250               CALL lbc_lnk( zitx1, 'T', -1. )   ;   CALL lbc_lnk( zity1, 'T', -1. )
2251            CASE( 'mixed oce-ice'        )
2252               SELECT CASE ( cp_ice_msh )
2253               CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
2254                  DO jj = 2, jpjm1
2255                     DO ji = fs_2, fs_jpim1   ! vector opt.
2256                        zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)   &
2257                           &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
2258                        zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)   &
2259                           &         + 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
2260                     END DO
2261                  END DO
2262               CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
2263                  DO jj = 2, jpjm1
2264                     DO ji = 2, jpim1   ! NO vector opt.
2265                        zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
2266                           &         + 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &
2267                           &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
2268                        zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
2269                           &         + 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &
2270                           &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
2271                     END DO
2272                  END DO
2273               CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
2274                  DO jj = 2, jpjm1
2275                     DO ji = 2, jpim1   ! NO vector opt.
2276                        zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
2277                           &         + 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &
2278                           &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
2279                        zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
2280                           &         + 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &
2281                           &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
2282                     END DO
2283                  END DO
2284               END SELECT
2285            END SELECT
2286            CALL lbc_lnk( zotx1, ssnd(jps_ocx1)%clgrid, -1. )   ;   CALL lbc_lnk( zoty1, ssnd(jps_ocy1)%clgrid, -1. )
2287            !
2288         ENDIF
2289         !
2290         !
2291         IF( TRIM( sn_snd_crt%clvor ) == 'eastward-northward' ) THEN             ! Rotation of the components
2292            !                                                                     ! Ocean component
2293            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->e', ztmp1 )       ! 1st component
2294            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->n', ztmp2 )       ! 2nd component
2295            zotx1(:,:) = ztmp1(:,:)                                                   ! overwrite the components
2296            zoty1(:,:) = ztmp2(:,:)
2297            IF( ssnd(jps_ivx1)%laction ) THEN                                     ! Ice component
2298               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 )    ! 1st component
2299               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 )    ! 2nd component
2300               zitx1(:,:) = ztmp1(:,:)                                                ! overwrite the components
2301               zity1(:,:) = ztmp2(:,:)
2302            ENDIF
2303         ENDIF
2304         !
2305         ! spherical coordinates to cartesian -> 2 components to 3 components
2306         IF( TRIM( sn_snd_crt%clvref ) == 'cartesian' ) THEN
2307            ztmp1(:,:) = zotx1(:,:)                     ! ocean currents
2308            ztmp2(:,:) = zoty1(:,:)
2309            CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 )
2310            !
2311            IF( ssnd(jps_ivx1)%laction ) THEN           ! ice velocities
2312               ztmp1(:,:) = zitx1(:,:)
2313               ztmp1(:,:) = zity1(:,:)
2314               CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 )
2315            ENDIF
2316         ENDIF
2317         !
2318         IF( ssnd(jps_ocx1)%laction )   CALL cpl_snd( jps_ocx1, isec, RESHAPE ( zotx1, (/jpi,jpj,1/) ), info )   ! ocean x current 1st grid
2319         IF( ssnd(jps_ocy1)%laction )   CALL cpl_snd( jps_ocy1, isec, RESHAPE ( zoty1, (/jpi,jpj,1/) ), info )   ! ocean y current 1st grid
2320         IF( ssnd(jps_ocz1)%laction )   CALL cpl_snd( jps_ocz1, isec, RESHAPE ( zotz1, (/jpi,jpj,1/) ), info )   ! ocean z current 1st grid
2321         !
2322         IF( ssnd(jps_ivx1)%laction )   CALL cpl_snd( jps_ivx1, isec, RESHAPE ( zitx1, (/jpi,jpj,1/) ), info )   ! ice   x current 1st grid
2323         IF( ssnd(jps_ivy1)%laction )   CALL cpl_snd( jps_ivy1, isec, RESHAPE ( zity1, (/jpi,jpj,1/) ), info )   ! ice   y current 1st grid
2324         IF( ssnd(jps_ivz1)%laction )   CALL cpl_snd( jps_ivz1, isec, RESHAPE ( zitz1, (/jpi,jpj,1/) ), info )   ! ice   z current 1st grid
2325         !
2326      ENDIF
2327      !
2328      !                                                      ! ------------------------- !
2329      !                                                      !  Surface current to waves !
2330      !                                                      ! ------------------------- !
2331      IF( ssnd(jps_ocxw)%laction .OR. ssnd(jps_ocyw)%laction ) THEN 
2332          !     
2333          !                                                  j+1  j     -----V---F
2334          ! surface velocity always sent from T point                    !       |
2335          !                                                       j      |   T   U
2336          !                                                              |       |
2337          !                                                   j   j-1   -I-------|
2338          !                                               (for I)        |       |
2339          !                                                             i-1  i   i
2340          !                                                              i      i+1 (for I)
2341          SELECT CASE( TRIM( sn_snd_crtw%cldes ) ) 
2342          CASE( 'oce only'             )      ! C-grid ==> T
2343             DO jj = 2, jpjm1 
2344                DO ji = fs_2, fs_jpim1   ! vector opt.
2345                   zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj  ,1) ) 
2346                   zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji , jj-1,1) ) 
2347                END DO
2348             END DO
2349          CASE( 'weighted oce and ice' )   
2350             SELECT CASE ( cp_ice_msh ) 
2351             CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
2352                DO jj = 2, jpjm1 
2353                   DO ji = fs_2, fs_jpim1   ! vector opt.
2354                      zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)   
2355                      zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj) 
2356                      zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj) 
2357                      zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj) 
2358                   END DO
2359                END DO
2360             CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
2361                DO jj = 2, jpjm1 
2362                   DO ji = 2, jpim1   ! NO vector opt.
2363                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   
2364                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   
2365                      zitx1(ji,jj) = 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     & 
2366                         &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj) 
2367                      zity1(ji,jj) = 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     & 
2368                         &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj) 
2369                   END DO
2370                END DO
2371             CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
2372                DO jj = 2, jpjm1 
2373                   DO ji = 2, jpim1   ! NO vector opt.
2374                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   
2375                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   
2376                      zitx1(ji,jj) = 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     & 
2377                         &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj) 
2378                      zity1(ji,jj) = 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     & 
2379                         &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj) 
2380                   END DO
2381                END DO
2382             END SELECT
2383             CALL lbc_lnk( zitx1, 'T', -1. )   ;   CALL lbc_lnk( zity1, 'T', -1. ) 
2384          CASE( 'mixed oce-ice'        ) 
2385             SELECT CASE ( cp_ice_msh ) 
2386             CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
2387                DO jj = 2, jpjm1 
2388                   DO ji = fs_2, fs_jpim1   ! vector opt.
2389                      zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)   & 
2390                         &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj) 
2391                      zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
2392                         &         + 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj) 
2393                   END DO
2394                END DO
2395             CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
2396                DO jj = 2, jpjm1 
2397                   DO ji = 2, jpim1   ! NO vector opt.
2398                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
2399                         &         + 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     & 
2400                         &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj) 
2401                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
2402                         &         + 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     & 
2403                         &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj) 
2404                   END DO
2405                END DO
2406             CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
2407                DO jj = 2, jpjm1 
2408                   DO ji = 2, jpim1   ! NO vector opt.
2409                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
2410                         &         + 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     & 
2411                         &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj) 
2412                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
2413                         &         + 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     & 
2414                         &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj) 
2415                   END DO
2416                END DO
2417             END SELECT
2418          END SELECT
2419         CALL lbc_lnk( zotx1, ssnd(jps_ocxw)%clgrid, -1. )   ;   CALL lbc_lnk( zoty1, ssnd(jps_ocyw)%clgrid, -1. ) 
2420         !
2421         !
2422         IF( TRIM( sn_snd_crtw%clvor ) == 'eastward-northward' ) THEN             ! Rotation of the components
2423         !                                                                        ! Ocean component
2424            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocxw)%clgrid, 'ij->e', ztmp1 )       ! 1st component 
2425            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocxw)%clgrid, 'ij->n', ztmp2 )       ! 2nd component 
2426            zotx1(:,:) = ztmp1(:,:)                                                   ! overwrite the components 
2427            zoty1(:,:) = ztmp2(:,:) 
2428            IF( ssnd(jps_ivx1)%laction ) THEN                                     ! Ice component
2429               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 )    ! 1st component 
2430               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 )    ! 2nd component 
2431               zitx1(:,:) = ztmp1(:,:)                                                ! overwrite the components 
2432               zity1(:,:) = ztmp2(:,:) 
2433            ENDIF
2434         ENDIF 
2435         !
2436!         ! spherical coordinates to cartesian -> 2 components to 3 components
2437!         IF( TRIM( sn_snd_crtw%clvref ) == 'cartesian' ) THEN
2438!            ztmp1(:,:) = zotx1(:,:)                     ! ocean currents
2439!            ztmp2(:,:) = zoty1(:,:)
2440!            CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 )
2441!            !
2442!            IF( ssnd(jps_ivx1)%laction ) THEN           ! ice velocities
2443!               ztmp1(:,:) = zitx1(:,:)
2444!               ztmp1(:,:) = zity1(:,:)
2445!               CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 )
2446!            ENDIF
2447!         ENDIF
2448         !
2449         IF( ssnd(jps_ocxw)%laction )   CALL cpl_snd( jps_ocxw, isec, RESHAPE ( zotx1, (/jpi,jpj,1/) ), info )   ! ocean x current 1st grid
2450         IF( ssnd(jps_ocyw)%laction )   CALL cpl_snd( jps_ocyw, isec, RESHAPE ( zoty1, (/jpi,jpj,1/) ), info )   ! ocean y current 1st grid
2451         
2452      ENDIF 
2453      !
2454      IF( ssnd(jps_ficet)%laction ) THEN
2455         CALL cpl_snd( jps_ficet, isec, RESHAPE ( fr_i, (/jpi,jpj,1/) ), info ) 
2456      END IF 
2457      !                                                      ! ------------------------- !
2458      !                                                      !   Water levels to waves   !
2459      !                                                      ! ------------------------- !
2460      IF( ssnd(jps_wlev)%laction ) THEN
2461         IF( ln_apr_dyn ) THEN 
2462            IF( kt /= nit000 ) THEN 
2463               ztmp1(:,:) = sshb(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) 
2464            ELSE 
2465               ztmp1(:,:) = sshb(:,:) 
2466            ENDIF 
2467         ELSE 
2468            ztmp1(:,:) = sshn(:,:) 
2469         ENDIF 
2470         CALL cpl_snd( jps_wlev  , isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info ) 
2471      END IF 
2472      !
2473      !  Fields sent by OPA to SAS when doing OPA<->SAS coupling
2474      !                                                        ! SSH
2475      IF( ssnd(jps_ssh )%laction )  THEN
2476         !                          ! removed inverse barometer ssh when Patm
2477         !                          forcing is used (for sea-ice dynamics)
2478         IF( ln_apr_dyn ) THEN   ;   ztmp1(:,:) = sshb(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
2479         ELSE                    ;   ztmp1(:,:) = sshn(:,:)
2480         ENDIF
2481         CALL cpl_snd( jps_ssh   , isec, RESHAPE ( ztmp1            , (/jpi,jpj,1/) ), info )
2482
2483      ENDIF
2484      !                                                        ! SSS
2485      IF( ssnd(jps_soce  )%laction )  THEN
2486         CALL cpl_snd( jps_soce  , isec, RESHAPE ( tsn(:,:,1,jp_sal), (/jpi,jpj,1/) ), info )
2487      ENDIF
2488      !                                                        ! first T level thickness
2489      IF( ssnd(jps_e3t1st )%laction )  THEN
2490         CALL cpl_snd( jps_e3t1st, isec, RESHAPE ( e3t_n(:,:,1)   , (/jpi,jpj,1/) ), info )
2491      ENDIF
2492      !                                                        ! Qsr fraction
2493      IF( ssnd(jps_fraqsr)%laction )  THEN
2494         CALL cpl_snd( jps_fraqsr, isec, RESHAPE ( fraqsr_1lev(:,:) , (/jpi,jpj,1/) ), info )
2495      ENDIF
2496      !
2497      !  Fields sent by SAS to OPA when OASIS coupling
2498      !                                                        ! Solar heat flux
2499      IF( ssnd(jps_qsroce)%laction )  CALL cpl_snd( jps_qsroce, isec, RESHAPE ( qsr , (/jpi,jpj,1/) ), info )
2500      IF( ssnd(jps_qnsoce)%laction )  CALL cpl_snd( jps_qnsoce, isec, RESHAPE ( qns , (/jpi,jpj,1/) ), info )
2501      IF( ssnd(jps_oemp  )%laction )  CALL cpl_snd( jps_oemp  , isec, RESHAPE ( emp , (/jpi,jpj,1/) ), info )
2502      IF( ssnd(jps_sflx  )%laction )  CALL cpl_snd( jps_sflx  , isec, RESHAPE ( sfx , (/jpi,jpj,1/) ), info )
2503      IF( ssnd(jps_otx1  )%laction )  CALL cpl_snd( jps_otx1  , isec, RESHAPE ( utau, (/jpi,jpj,1/) ), info )
2504      IF( ssnd(jps_oty1  )%laction )  CALL cpl_snd( jps_oty1  , isec, RESHAPE ( vtau, (/jpi,jpj,1/) ), info )
2505      IF( ssnd(jps_rnf   )%laction )  CALL cpl_snd( jps_rnf   , isec, RESHAPE ( rnf , (/jpi,jpj,1/) ), info )
2506      IF( ssnd(jps_taum  )%laction )  CALL cpl_snd( jps_taum  , isec, RESHAPE ( taum, (/jpi,jpj,1/) ), info )
2507
2508      CALL wrk_dealloc( jpi,jpj,       zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
2509      CALL wrk_dealloc( jpi,jpj,jpl,   ztmp3, ztmp4 )
2510      !
2511      IF( nn_timing == 1 )   CALL timing_stop('sbc_cpl_snd')
2512      !
2513   END SUBROUTINE sbc_cpl_snd
2514   
2515   !!======================================================================
2516END MODULE sbccpl
Note: See TracBrowser for help on using the repository browser.