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

source: branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 8752

Last change on this file since 8752 was 8752, checked in by dancopsey, 6 years ago

Merged in main ICEMODEL branch (branches/2017/dev_r8183_ICEMODEL) from revision 8587 to 8726.

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