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

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

source: branches/UKMO/dev_r5518_GO6_fix_zemp_ice/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 11055

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

Put the switch back in again.

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