source: branches/UKMO/dev_r5518_cleanup_1d_cpl/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 10046

Last change on this file since 10046 was 10046, checked in by dancopsey, 2 years ago

Re-enable 2D coupling of icesheet mass. To switch between 2D to 0D coupling of ice mass change sn_rcv_grnm%cldes from 'coupled' to 'coupled0d' (with a similar change to sn_rcv_antm).

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