source: branches/UKMO/dev_r5518_GSI7_GSI8_landice_bitcomp_medusa/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 6659

Last change on this file since 6659 was 6659, checked in by frrh, 4 years ago

Merge branch dev_r5518_coupling_GSI7_GSI8_landice_bitcomp from
revision 6363 to 6651 inclusive.

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