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

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

source: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 7990

Last change on this file since 7990 was 7990, checked in by gm, 7 years ago

#1880 (HPC-09): OPA remove avmu, avmv from zdf modules + move CALL tke(gls)_rst & gls_rst in zdftke(gls) + rename zdftmx and zdfqiao

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