New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sbccpl.F90 in branches/UKMO/r8727_WAVE-2_Clementi_add_coupling/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

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

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

First set of changes for ticket #1980
Addition of the Phillips parametrization for vertical Stokes drift

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