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_cleanup_1d_cpl/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

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

Last change on this file since 9280 was 9280, checked in by frrh, 6 years ago

Add support for reading in nn_cpl_river, the number of rivers to be dealt with
in atmos-ocean coupling.

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