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

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

source: branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 12434

Last change on this file since 12434 was 12434, checked in by jcastill, 4 years ago

Merge development branch AMM15_v3_6_STABLE_package_collate_utils325, see Met Office utils branch 325. As per reviewer's comments, indentation has been included for the added if sentence

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