source: branches/UKMO/dev_r7750_GO6_package_oasis_timers/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 8887

Last change on this file since 8887 was 8887, checked in by andmirek, 3 years ago

#1978 coupling timers

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