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

source: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 8813

Last change on this file since 8813 was 8813, checked in by gm, 6 years ago

#1911 (ENHANCE-09): PART I.3 - phasing with updated branch dev_r8183_ICEMODEL revision 8787

  • 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      !                                                      ! ========================= !
1946      !                                                      !      Transmitted Qsr      !   [W/m2]
1947      !                                                      ! ========================= !
1948      SELECT CASE( nice_jules )
1949      CASE( np_jules_OFF    )       !==  No Jules coupler  ==!
1950         !
1951!!gm         ! former coding was   
1952!!gm         ! Coupled case: since cloud cover is not received from atmosphere
1953!!gm         !               ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
1954!!gm         !     fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
1955!!gm         !     fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
1956!!gm         
1957!!gm         ! to retrieve that coding, we needed to access h_i & h_s from here
1958!!gm         ! we could even retrieve cloud fraction from the coupler
1959!!gm         !
1960!!gm         zfrqsr_tr_i(:,:,:) = 0._wp   !   surface transmission parameter
1961!!gm         !
1962!!gm         DO jl = 1, jpl
1963!!gm            DO jj = 1 , jpj
1964!!gm               DO ji = 1, jpi
1965!!gm                  !              !--- surface transmission parameter (Grenfell Maykut 77) --- !
1966!!gm                  zfrqsr_tr_i(ji,jj,jl) = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice
1967!!gm                  !
1968!!gm                  !              ! --- influence of snow and thin ice --- !
1969!!gm                  IF ( phs(ji,jj,jl) >= 0.0_wp )   zfrqsr_tr_i(ji,jj,jl) = 0._wp   !   snow fully opaque
1970!!gm                  IF ( phi(ji,jj,jl) <= 0.1_wp )   zfrqsr_tr_i(ji,jj,jl) = 1._wp   !   thin ice transmits all solar radiation
1971!!gm               END DO
1972!!gm            END DO
1973!!gm         END DO
1974!!gm         !
1975!!gm         qsr_ice_tr(:,:,:) =   zfrqsr_tr_i(:,:,:) * qsr_ice(:,:,:)               !   transmitted solar radiation
1976!!gm         !
1977!!gm better coding of the above calculation:
1978         !
1979         !                    ! ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
1980         ztri = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice    ! surface transmission parameter (Grenfell Maykut 77)
1981         !
1982         qsr_ice_tr(:,:,:) = ztri * qsr_ice(:,:,:)
1983         WHERE( phs(:,:,:) >= 0.0_wp )   qsr_ice_tr(:,:,:) = 0._wp            ! snow fully opaque
1984         WHERE( phi(:,:,:) <= 0.1_wp )   qsr_ice_tr(:,:,:) = qsr_ice(:,:,:)   ! thin ice transmits all solar radiation
1985!!gm end
1986         !     
1987      CASE( np_jules_ACTIVE )       !==  Jules coupler is active  ==!
1988         !
1989         !                    ! ===> here we must receive the qsr_ice_tr array from the coupler
1990         !                           for now just assume zero (fully opaque ice)
1991         qsr_ice_tr(:,:,:) = 0._wp
1992         !
1993      END SELECT
1994      !
1995      IF( nn_timing == 1 )   CALL timing_stop('sbc_cpl_ice_flx')
1996      !
1997   END SUBROUTINE sbc_cpl_ice_flx
1998   
1999   
2000   SUBROUTINE sbc_cpl_snd( kt )
2001      !!----------------------------------------------------------------------
2002      !!             ***  ROUTINE sbc_cpl_snd  ***
2003      !!
2004      !! ** Purpose :   provide the ocean-ice informations to the atmosphere
2005      !!
2006      !! ** Method  :   send to the atmosphere through a call to cpl_snd
2007      !!              all the needed fields (as defined in sbc_cpl_init)
2008      !!----------------------------------------------------------------------
2009      INTEGER, INTENT(in) ::   kt
2010      !
2011      INTEGER ::   ji, jj, jl   ! dummy loop indices
2012      INTEGER ::   isec, info   ! local integer
2013      REAL(wp) ::   zumax, zvmax
2014      REAL(wp), POINTER, DIMENSION(:,:)   ::   zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1
2015      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmp3, ztmp4   
2016      !!----------------------------------------------------------------------
2017      !
2018      IF( nn_timing == 1 )   CALL timing_start('sbc_cpl_snd')
2019      !
2020      CALL wrk_alloc( jpi,jpj,       zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
2021      CALL wrk_alloc( jpi,jpj,jpl,   ztmp3, ztmp4 )
2022
2023      isec = ( kt - nit000 ) * NINT( rdt )        ! date of exchanges
2024
2025      zfr_l(:,:) = 1.- fr_i(:,:)
2026      !                                                      ! ------------------------- !
2027      !                                                      !    Surface temperature    !   in Kelvin
2028      !                                                      ! ------------------------- !
2029      IF( ssnd(jps_toce)%laction .OR. ssnd(jps_tice)%laction .OR. ssnd(jps_tmix)%laction ) THEN
2030         
2031         IF ( nn_components == jp_iam_opa ) THEN
2032            ztmp1(:,:) = tsn(:,:,1,jp_tem)   ! send temperature as it is (potential or conservative) -> use of l_useCT on the received part
2033         ELSE
2034            ! we must send the surface potential temperature
2035            IF( l_useCT )  THEN    ;   ztmp1(:,:) = eos_pt_from_ct( tsn(:,:,1,jp_tem), tsn(:,:,1,jp_sal) )
2036            ELSE                   ;   ztmp1(:,:) = tsn(:,:,1,jp_tem)
2037            ENDIF
2038            !
2039            SELECT CASE( sn_snd_temp%cldes)
2040            CASE( 'oce only'             )   ;   ztmp1(:,:) =   ztmp1(:,:) + rt0
2041            CASE( 'oce and ice'          )   ;   ztmp1(:,:) =   ztmp1(:,:) + rt0
2042               SELECT CASE( sn_snd_temp%clcat )
2043               CASE( 'yes' )   
2044                  ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl)
2045               CASE( 'no' )
2046                  WHERE( SUM( a_i, dim=3 ) /= 0. )
2047                     ztmp3(:,:,1) = SUM( tn_ice * a_i, dim=3 ) / SUM( a_i, dim=3 )
2048                  ELSEWHERE
2049                     ztmp3(:,:,1) = rt0
2050                  END WHERE
2051               CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
2052               END SELECT
2053            CASE( 'weighted oce and ice' )   ;   ztmp1(:,:) = ( ztmp1(:,:) + rt0 ) * zfr_l(:,:)   
2054               SELECT CASE( sn_snd_temp%clcat )
2055               CASE( 'yes' )   
2056                  ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
2057               CASE( 'no' )
2058                  ztmp3(:,:,:) = 0.0
2059                  DO jl=1,jpl
2060                     ztmp3(:,:,1) = ztmp3(:,:,1) + tn_ice(:,:,jl) * a_i(:,:,jl)
2061                  ENDDO
2062               CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
2063               END SELECT
2064            CASE( 'mixed oce-ice'        )   
2065               ztmp1(:,:) = ( ztmp1(:,:) + rt0 ) * zfr_l(:,:) 
2066               DO jl=1,jpl
2067                  ztmp1(:,:) = ztmp1(:,:) + tn_ice(:,:,jl) * a_i(:,:,jl)
2068               ENDDO
2069            CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%cldes' )
2070            END SELECT
2071         ENDIF
2072         IF( ssnd(jps_toce)%laction )   CALL cpl_snd( jps_toce, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
2073         IF( ssnd(jps_tice)%laction )   CALL cpl_snd( jps_tice, isec, ztmp3, info )
2074         IF( ssnd(jps_tmix)%laction )   CALL cpl_snd( jps_tmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
2075      ENDIF
2076      !                                                      ! ------------------------- !
2077      !                                                      !           Albedo          !
2078      !                                                      ! ------------------------- !
2079      IF( ssnd(jps_albice)%laction ) THEN                         ! ice
2080          SELECT CASE( sn_snd_alb%cldes )
2081          CASE( 'ice' )
2082             SELECT CASE( sn_snd_alb%clcat )
2083             CASE( 'yes' )   
2084                ztmp3(:,:,1:jpl) = alb_ice(:,:,1:jpl)
2085             CASE( 'no' )
2086                WHERE( SUM( a_i, dim=3 ) /= 0. )
2087                   ztmp1(:,:) = SUM( alb_ice (:,:,1:jpl) * a_i(:,:,1:jpl), dim=3 ) / SUM( a_i(:,:,1:jpl), dim=3 )
2088                ELSEWHERE
2089                   ztmp1(:,:) = alb_oce_mix(:,:)
2090                END WHERE
2091             CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_alb%clcat' )
2092             END SELECT
2093          CASE( 'weighted ice' )   ;
2094             SELECT CASE( sn_snd_alb%clcat )
2095             CASE( 'yes' )   
2096                ztmp3(:,:,1:jpl) =  alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
2097             CASE( 'no' )
2098                WHERE( fr_i (:,:) > 0. )
2099                   ztmp1(:,:) = SUM (  alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl), dim=3 )
2100                ELSEWHERE
2101                   ztmp1(:,:) = 0.
2102                END WHERE
2103             CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_ice%clcat' )
2104             END SELECT
2105          CASE default      ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_alb%cldes' )
2106         END SELECT
2107
2108         SELECT CASE( sn_snd_alb%clcat )
2109            CASE( 'yes' )   
2110               CALL cpl_snd( jps_albice, isec, ztmp3, info )      !-> MV this has never been checked in coupled mode
2111            CASE( 'no'  )   
2112               CALL cpl_snd( jps_albice, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info ) 
2113         END SELECT
2114      ENDIF
2115
2116      IF( ssnd(jps_albmix)%laction ) THEN                         ! mixed ice-ocean
2117         ztmp1(:,:) = alb_oce_mix(:,:) * zfr_l(:,:)
2118         DO jl = 1, jpl
2119            ztmp1(:,:) = ztmp1(:,:) + alb_ice(:,:,jl) * a_i(:,:,jl)
2120         END DO
2121         CALL cpl_snd( jps_albmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
2122      ENDIF
2123      !                                                      ! ------------------------- !
2124      !                                                      !  Ice fraction & Thickness !
2125      !                                                      ! ------------------------- !
2126      ! Send ice fraction field to atmosphere
2127      IF( ssnd(jps_fice)%laction ) THEN
2128         SELECT CASE( sn_snd_thick%clcat )
2129         CASE( 'yes' )   ;   ztmp3(:,:,1:jpl) =  a_i(:,:,1:jpl)
2130         CASE( 'no'  )   ;   ztmp3(:,:,1    ) = fr_i(:,:      )
2131         CASE default    ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
2132         END SELECT
2133         IF( ssnd(jps_fice)%laction )   CALL cpl_snd( jps_fice, isec, ztmp3, info )
2134      ENDIF
2135     
2136      ! Send ice fraction field to OPA (sent by SAS in SAS-OPA coupling)
2137      IF( ssnd(jps_fice2)%laction ) THEN
2138         ztmp3(:,:,1) = fr_i(:,:)
2139         IF( ssnd(jps_fice2)%laction )   CALL cpl_snd( jps_fice2, isec, ztmp3, info )
2140      ENDIF
2141
2142      ! Send ice and snow thickness field
2143      IF( ssnd(jps_hice)%laction .OR. ssnd(jps_hsnw)%laction ) THEN
2144         SELECT CASE( sn_snd_thick%cldes)
2145         CASE( 'none'                  )       ! nothing to do
2146         CASE( 'weighted ice and snow' )   
2147            SELECT CASE( sn_snd_thick%clcat )
2148            CASE( 'yes' )   
2149               ztmp3(:,:,1:jpl) =  h_i(:,:,1:jpl) * a_i(:,:,1:jpl)
2150               ztmp4(:,:,1:jpl) =  h_s(:,:,1:jpl) * a_i(:,:,1:jpl)
2151            CASE( 'no' )
2152               ztmp3(:,:,:) = 0.0   ;  ztmp4(:,:,:) = 0.0
2153               DO jl=1,jpl
2154                  ztmp3(:,:,1) = ztmp3(:,:,1) + h_i(:,:,jl) * a_i(:,:,jl)
2155                  ztmp4(:,:,1) = ztmp4(:,:,1) + h_s(:,:,jl) * a_i(:,:,jl)
2156               ENDDO
2157            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
2158            END SELECT
2159         CASE( 'ice and snow'         )   
2160            SELECT CASE( sn_snd_thick%clcat )
2161            CASE( 'yes' )
2162               ztmp3(:,:,1:jpl) = h_i(:,:,1:jpl)
2163               ztmp4(:,:,1:jpl) = h_s(:,:,1:jpl)
2164            CASE( 'no' )
2165               WHERE( SUM( a_i, dim=3 ) /= 0. )
2166                  ztmp3(:,:,1) = SUM( h_i * a_i, dim=3 ) / SUM( a_i, dim=3 )
2167                  ztmp4(:,:,1) = SUM( h_s * a_i, dim=3 ) / SUM( a_i, dim=3 )
2168               ELSEWHERE
2169                 ztmp3(:,:,1) = 0.
2170                 ztmp4(:,:,1) = 0.
2171               END WHERE
2172            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
2173            END SELECT
2174         CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%cldes' )
2175         END SELECT
2176         IF( ssnd(jps_hice)%laction )   CALL cpl_snd( jps_hice, isec, ztmp3, info )
2177         IF( ssnd(jps_hsnw)%laction )   CALL cpl_snd( jps_hsnw, isec, ztmp4, info )
2178      ENDIF
2179      !                                                      ! ------------------------- !
2180      !                                                      !  CO2 flux from PISCES     !
2181      !                                                      ! ------------------------- !
2182      IF( ssnd(jps_co2)%laction .AND. l_co2cpl )   CALL cpl_snd( jps_co2, isec, RESHAPE ( oce_co2, (/jpi,jpj,1/) ) , info )
2183      !
2184      !                                                      ! ------------------------- !
2185      IF( ssnd(jps_ocx1)%laction ) THEN                      !      Surface current      !
2186         !                                                   ! ------------------------- !
2187         !   
2188         !                                                  j+1   j     -----V---F
2189         ! surface velocity always sent from T point                     !       |
2190         !                                                        j      |   T   U
2191         !                                                               |       |
2192         !                                                   j    j-1   -I-------|
2193         !                                               (for I)         |       |
2194         !                                                              i-1  i   i
2195         !                                                               i      i+1 (for I)
2196         IF( nn_components == jp_iam_opa ) THEN
2197            zotx1(:,:) = un(:,:,1) 
2198            zoty1(:,:) = vn(:,:,1) 
2199         ELSE       
2200            SELECT CASE( TRIM( sn_snd_crt%cldes ) )
2201            CASE( 'oce only'             )      ! C-grid ==> T
2202               DO jj = 2, jpjm1
2203                  DO ji = fs_2, fs_jpim1   ! vector opt.
2204                     zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj  ,1) )
2205                     zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji  ,jj-1,1) ) 
2206                  END DO
2207               END DO
2208            CASE( 'weighted oce and ice' )   
2209               SELECT CASE ( cp_ice_msh )
2210               CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
2211                  DO jj = 2, jpjm1
2212                     DO ji = fs_2, fs_jpim1   ! vector opt.
2213                        zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj) 
2214                        zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)
2215                        zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
2216                        zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
2217                     END DO
2218                  END DO
2219               CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
2220                  DO jj = 2, jpjm1
2221                     DO ji = 2, jpim1   ! NO vector opt.
2222                        zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj) 
2223                        zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj) 
2224                        zitx1(ji,jj) = 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &
2225                           &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
2226                        zity1(ji,jj) = 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &
2227                           &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
2228                     END DO
2229                  END DO
2230               CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
2231                  DO jj = 2, jpjm1
2232                     DO ji = 2, jpim1   ! NO vector opt.
2233                        zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj) 
2234                        zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj) 
2235                        zitx1(ji,jj) = 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &
2236                           &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
2237                        zity1(ji,jj) = 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &
2238                           &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
2239                     END DO
2240                  END DO
2241               END SELECT
2242               CALL lbc_lnk( zitx1, 'T', -1. )   ;   CALL lbc_lnk( zity1, 'T', -1. )
2243            CASE( 'mixed oce-ice'        )
2244               SELECT CASE ( cp_ice_msh )
2245               CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
2246                  DO jj = 2, jpjm1
2247                     DO ji = fs_2, fs_jpim1   ! vector opt.
2248                        zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)   &
2249                           &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
2250                        zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)   &
2251                           &         + 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
2252                     END DO
2253                  END DO
2254               CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
2255                  DO jj = 2, jpjm1
2256                     DO ji = 2, jpim1   ! NO vector opt.
2257                        zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
2258                           &         + 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &
2259                           &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
2260                        zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
2261                           &         + 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &
2262                           &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
2263                     END DO
2264                  END DO
2265               CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
2266                  DO jj = 2, jpjm1
2267                     DO ji = 2, jpim1   ! NO vector opt.
2268                        zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
2269                           &         + 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &
2270                           &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
2271                        zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
2272                           &         + 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &
2273                           &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
2274                     END DO
2275                  END DO
2276               END SELECT
2277            END SELECT
2278            CALL lbc_lnk( zotx1, ssnd(jps_ocx1)%clgrid, -1. )   ;   CALL lbc_lnk( zoty1, ssnd(jps_ocy1)%clgrid, -1. )
2279            !
2280         ENDIF
2281         !
2282         !
2283         IF( TRIM( sn_snd_crt%clvor ) == 'eastward-northward' ) THEN             ! Rotation of the components
2284            !                                                                     ! Ocean component
2285            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->e', ztmp1 )       ! 1st component
2286            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->n', ztmp2 )       ! 2nd component
2287            zotx1(:,:) = ztmp1(:,:)                                                   ! overwrite the components
2288            zoty1(:,:) = ztmp2(:,:)
2289            IF( ssnd(jps_ivx1)%laction ) THEN                                     ! Ice component
2290               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 )    ! 1st component
2291               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 )    ! 2nd component
2292               zitx1(:,:) = ztmp1(:,:)                                                ! overwrite the components
2293               zity1(:,:) = ztmp2(:,:)
2294            ENDIF
2295         ENDIF
2296         !
2297         ! spherical coordinates to cartesian -> 2 components to 3 components
2298         IF( TRIM( sn_snd_crt%clvref ) == 'cartesian' ) THEN
2299            ztmp1(:,:) = zotx1(:,:)                     ! ocean currents
2300            ztmp2(:,:) = zoty1(:,:)
2301            CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 )
2302            !
2303            IF( ssnd(jps_ivx1)%laction ) THEN           ! ice velocities
2304               ztmp1(:,:) = zitx1(:,:)
2305               ztmp1(:,:) = zity1(:,:)
2306               CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 )
2307            ENDIF
2308         ENDIF
2309         !
2310         IF( ssnd(jps_ocx1)%laction )   CALL cpl_snd( jps_ocx1, isec, RESHAPE ( zotx1, (/jpi,jpj,1/) ), info )   ! ocean x current 1st grid
2311         IF( ssnd(jps_ocy1)%laction )   CALL cpl_snd( jps_ocy1, isec, RESHAPE ( zoty1, (/jpi,jpj,1/) ), info )   ! ocean y current 1st grid
2312         IF( ssnd(jps_ocz1)%laction )   CALL cpl_snd( jps_ocz1, isec, RESHAPE ( zotz1, (/jpi,jpj,1/) ), info )   ! ocean z current 1st grid
2313         !
2314         IF( ssnd(jps_ivx1)%laction )   CALL cpl_snd( jps_ivx1, isec, RESHAPE ( zitx1, (/jpi,jpj,1/) ), info )   ! ice   x current 1st grid
2315         IF( ssnd(jps_ivy1)%laction )   CALL cpl_snd( jps_ivy1, isec, RESHAPE ( zity1, (/jpi,jpj,1/) ), info )   ! ice   y current 1st grid
2316         IF( ssnd(jps_ivz1)%laction )   CALL cpl_snd( jps_ivz1, isec, RESHAPE ( zitz1, (/jpi,jpj,1/) ), info )   ! ice   z current 1st grid
2317         !
2318      ENDIF
2319      !
2320      !                                                      ! ------------------------- !
2321      !                                                      !  Surface current to waves !
2322      !                                                      ! ------------------------- !
2323      IF( ssnd(jps_ocxw)%laction .OR. ssnd(jps_ocyw)%laction ) THEN 
2324          !     
2325          !                                                  j+1  j     -----V---F
2326          ! surface velocity always sent from T point                    !       |
2327          !                                                       j      |   T   U
2328          !                                                              |       |
2329          !                                                   j   j-1   -I-------|
2330          !                                               (for I)        |       |
2331          !                                                             i-1  i   i
2332          !                                                              i      i+1 (for I)
2333          SELECT CASE( TRIM( sn_snd_crtw%cldes ) ) 
2334          CASE( 'oce only'             )      ! C-grid ==> T
2335             DO jj = 2, jpjm1 
2336                DO ji = fs_2, fs_jpim1   ! vector opt.
2337                   zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj  ,1) ) 
2338                   zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji , jj-1,1) ) 
2339                END DO
2340             END DO
2341          CASE( 'weighted oce and ice' )   
2342             SELECT CASE ( cp_ice_msh ) 
2343             CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
2344                DO jj = 2, jpjm1 
2345                   DO ji = fs_2, fs_jpim1   ! vector opt.
2346                      zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)   
2347                      zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj) 
2348                      zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj) 
2349                      zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj) 
2350                   END DO
2351                END DO
2352             CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
2353                DO jj = 2, jpjm1 
2354                   DO ji = 2, jpim1   ! NO vector opt.
2355                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   
2356                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   
2357                      zitx1(ji,jj) = 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     & 
2358                         &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj) 
2359                      zity1(ji,jj) = 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     & 
2360                         &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj) 
2361                   END DO
2362                END DO
2363             CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
2364                DO jj = 2, jpjm1 
2365                   DO ji = 2, jpim1   ! NO vector opt.
2366                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   
2367                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   
2368                      zitx1(ji,jj) = 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     & 
2369                         &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj) 
2370                      zity1(ji,jj) = 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     & 
2371                         &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj) 
2372                   END DO
2373                END DO
2374             END SELECT
2375             CALL lbc_lnk( zitx1, 'T', -1. )   ;   CALL lbc_lnk( zity1, 'T', -1. ) 
2376          CASE( 'mixed oce-ice'        ) 
2377             SELECT CASE ( cp_ice_msh ) 
2378             CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
2379                DO jj = 2, jpjm1 
2380                   DO ji = fs_2, fs_jpim1   ! vector opt.
2381                      zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)   & 
2382                         &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj) 
2383                      zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
2384                         &         + 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj) 
2385                   END DO
2386                END DO
2387             CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
2388                DO jj = 2, jpjm1 
2389                   DO ji = 2, jpim1   ! NO vector opt.
2390                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
2391                         &         + 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     & 
2392                         &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj) 
2393                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
2394                         &         + 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     & 
2395                         &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj) 
2396                   END DO
2397                END DO
2398             CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
2399                DO jj = 2, jpjm1 
2400                   DO ji = 2, jpim1   ! NO vector opt.
2401                      zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
2402                         &         + 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     & 
2403                         &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj) 
2404                      zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
2405                         &         + 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     & 
2406                         &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj) 
2407                   END DO
2408                END DO
2409             END SELECT
2410          END SELECT
2411         CALL lbc_lnk( zotx1, ssnd(jps_ocxw)%clgrid, -1. )   ;   CALL lbc_lnk( zoty1, ssnd(jps_ocyw)%clgrid, -1. ) 
2412         !
2413         !
2414         IF( TRIM( sn_snd_crtw%clvor ) == 'eastward-northward' ) THEN             ! Rotation of the components
2415         !                                                                        ! Ocean component
2416            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocxw)%clgrid, 'ij->e', ztmp1 )       ! 1st component 
2417            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocxw)%clgrid, 'ij->n', ztmp2 )       ! 2nd component 
2418            zotx1(:,:) = ztmp1(:,:)                                                   ! overwrite the components 
2419            zoty1(:,:) = ztmp2(:,:) 
2420            IF( ssnd(jps_ivx1)%laction ) THEN                                     ! Ice component
2421               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 )    ! 1st component 
2422               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 )    ! 2nd component 
2423               zitx1(:,:) = ztmp1(:,:)                                                ! overwrite the components 
2424               zity1(:,:) = ztmp2(:,:) 
2425            ENDIF
2426         ENDIF 
2427         !
2428!         ! spherical coordinates to cartesian -> 2 components to 3 components
2429!         IF( TRIM( sn_snd_crtw%clvref ) == 'cartesian' ) THEN
2430!            ztmp1(:,:) = zotx1(:,:)                     ! ocean currents
2431!            ztmp2(:,:) = zoty1(:,:)
2432!            CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 )
2433!            !
2434!            IF( ssnd(jps_ivx1)%laction ) THEN           ! ice velocities
2435!               ztmp1(:,:) = zitx1(:,:)
2436!               ztmp1(:,:) = zity1(:,:)
2437!               CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 )
2438!            ENDIF
2439!         ENDIF
2440         !
2441         IF( ssnd(jps_ocxw)%laction )   CALL cpl_snd( jps_ocxw, isec, RESHAPE ( zotx1, (/jpi,jpj,1/) ), info )   ! ocean x current 1st grid
2442         IF( ssnd(jps_ocyw)%laction )   CALL cpl_snd( jps_ocyw, isec, RESHAPE ( zoty1, (/jpi,jpj,1/) ), info )   ! ocean y current 1st grid
2443         
2444      ENDIF 
2445      !
2446      IF( ssnd(jps_ficet)%laction ) THEN
2447         CALL cpl_snd( jps_ficet, isec, RESHAPE ( fr_i, (/jpi,jpj,1/) ), info ) 
2448      END IF 
2449      !                                                      ! ------------------------- !
2450      !                                                      !   Water levels to waves   !
2451      !                                                      ! ------------------------- !
2452      IF( ssnd(jps_wlev)%laction ) THEN
2453         IF( ln_apr_dyn ) THEN 
2454            IF( kt /= nit000 ) THEN 
2455               ztmp1(:,:) = sshb(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) 
2456            ELSE 
2457               ztmp1(:,:) = sshb(:,:) 
2458            ENDIF 
2459         ELSE 
2460            ztmp1(:,:) = sshn(:,:) 
2461         ENDIF 
2462         CALL cpl_snd( jps_wlev  , isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info ) 
2463      END IF 
2464      !
2465      !  Fields sent by OPA to SAS when doing OPA<->SAS coupling
2466      !                                                        ! SSH
2467      IF( ssnd(jps_ssh )%laction )  THEN
2468         !                          ! removed inverse barometer ssh when Patm
2469         !                          forcing is used (for sea-ice dynamics)
2470         IF( ln_apr_dyn ) THEN   ;   ztmp1(:,:) = sshb(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
2471         ELSE                    ;   ztmp1(:,:) = sshn(:,:)
2472         ENDIF
2473         CALL cpl_snd( jps_ssh   , isec, RESHAPE ( ztmp1            , (/jpi,jpj,1/) ), info )
2474
2475      ENDIF
2476      !                                                        ! SSS
2477      IF( ssnd(jps_soce  )%laction )  THEN
2478         CALL cpl_snd( jps_soce  , isec, RESHAPE ( tsn(:,:,1,jp_sal), (/jpi,jpj,1/) ), info )
2479      ENDIF
2480      !                                                        ! first T level thickness
2481      IF( ssnd(jps_e3t1st )%laction )  THEN
2482         CALL cpl_snd( jps_e3t1st, isec, RESHAPE ( e3t_n(:,:,1)   , (/jpi,jpj,1/) ), info )
2483      ENDIF
2484      !                                                        ! Qsr fraction
2485      IF( ssnd(jps_fraqsr)%laction )  THEN
2486         CALL cpl_snd( jps_fraqsr, isec, RESHAPE ( fraqsr_1lev(:,:) , (/jpi,jpj,1/) ), info )
2487      ENDIF
2488      !
2489      !  Fields sent by SAS to OPA when OASIS coupling
2490      !                                                        ! Solar heat flux
2491      IF( ssnd(jps_qsroce)%laction )  CALL cpl_snd( jps_qsroce, isec, RESHAPE ( qsr , (/jpi,jpj,1/) ), info )
2492      IF( ssnd(jps_qnsoce)%laction )  CALL cpl_snd( jps_qnsoce, isec, RESHAPE ( qns , (/jpi,jpj,1/) ), info )
2493      IF( ssnd(jps_oemp  )%laction )  CALL cpl_snd( jps_oemp  , isec, RESHAPE ( emp , (/jpi,jpj,1/) ), info )
2494      IF( ssnd(jps_sflx  )%laction )  CALL cpl_snd( jps_sflx  , isec, RESHAPE ( sfx , (/jpi,jpj,1/) ), info )
2495      IF( ssnd(jps_otx1  )%laction )  CALL cpl_snd( jps_otx1  , isec, RESHAPE ( utau, (/jpi,jpj,1/) ), info )
2496      IF( ssnd(jps_oty1  )%laction )  CALL cpl_snd( jps_oty1  , isec, RESHAPE ( vtau, (/jpi,jpj,1/) ), info )
2497      IF( ssnd(jps_rnf   )%laction )  CALL cpl_snd( jps_rnf   , isec, RESHAPE ( rnf , (/jpi,jpj,1/) ), info )
2498      IF( ssnd(jps_taum  )%laction )  CALL cpl_snd( jps_taum  , isec, RESHAPE ( taum, (/jpi,jpj,1/) ), info )
2499
2500      CALL wrk_dealloc( jpi,jpj,       zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
2501      CALL wrk_dealloc( jpi,jpj,jpl,   ztmp3, ztmp4 )
2502      !
2503      IF( nn_timing == 1 )   CALL timing_stop('sbc_cpl_snd')
2504      !
2505   END SUBROUTINE sbc_cpl_snd
2506   
2507   !!======================================================================
2508END MODULE sbccpl
Note: See TracBrowser for help on using the repository browser.