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/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 8882

Last change on this file since 8882 was 8882, checked in by flavoni, 6 years ago

dev_CNRS_2017 branch: merged dev_r7881_ENHANCE09_RK3 with trunk r8864

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