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

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

Last change on this file since 9261 was 9261, checked in by andmirek, 6 years ago

#1978 only 1 and 2 for nn_timing

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