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

source: branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 8136

Last change on this file since 8136 was 8136, checked in by frrh, 7 years ago

Remove temporary updates for test purposes to enable live testing of chloro
coupling.

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