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

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

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 8726

Last change on this file since 8726 was 8726, checked in by vancop, 6 years ago

Jules coupler functionalities

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