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

source: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 5377

Last change on this file since 5377 was 5377, checked in by dancopsey, 9 years ago

Added sn_snd_mpnd to the list of namelist variables to load.

File size: 100.0 KB
Line 
1MODULE sbccpl
2   !!======================================================================
3   !!                       ***  MODULE  sbccpl  ***
4   !! Surface Boundary Condition :  momentum, heat and freshwater fluxes in coupled mode
5   !!======================================================================
6   !! History :  2.0  ! 2007-06  (R. Redler, N. Keenlyside, W. Park) Original code split into flxmod & taumod
7   !!            3.0  ! 2008-02  (G. Madec, C Talandier)  surface module
8   !!            3.1  ! 2009_02  (G. Madec, S. Masson, E. Maisonave, A. Caubel) generic coupled interface
9   !!            3.4  ! 2011_11  (C. Harris) more flexibility + multi-category fields
10   !!----------------------------------------------------------------------
11   !!----------------------------------------------------------------------
12   !!   namsbc_cpl      : coupled formulation namlist
13   !!   sbc_cpl_init    : initialisation of the coupled exchanges
14   !!   sbc_cpl_rcv     : receive fields from the atmosphere over the ocean (ocean only)
15   !!                     receive stress from the atmosphere over the ocean (ocean-ice case)
16   !!   sbc_cpl_ice_tau : receive stress from the atmosphere over ice
17   !!   sbc_cpl_ice_flx : receive fluxes from the atmosphere over ice
18   !!   sbc_cpl_snd     : send     fields to the atmosphere
19   !!----------------------------------------------------------------------
20   USE dom_oce         ! ocean space and time domain
21   USE sbc_oce         ! Surface boundary condition: ocean fields
22   USE sbc_ice         ! Surface boundary condition: ice fields
23   USE sbcdcy          ! surface boundary condition: diurnal cycle
24   USE phycst          ! physical constants
25#if defined key_lim3
26   USE par_ice         ! ice parameters
27   USE ice             ! ice variables
28#endif
29#if defined key_lim2
30   USE par_ice_2       ! ice parameters
31   USE ice_2           ! ice variables
32#endif
33   USE cpl_oasis3      ! OASIS3 coupling
34   USE geo2ocean       !
35   USE oce   , ONLY : tsn, un, vn
36   USE albedo          !
37   USE in_out_manager  ! I/O manager
38   USE iom             ! NetCDF library
39   USE lib_mpp         ! distribued memory computing library
40   USE wrk_nemo        ! work arrays
41   USE timing          ! Timing
42   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
43#if defined key_cpl_carbon_cycle
44   USE p4zflx, ONLY : oce_co2
45#endif
46   IMPLICIT NONE
47   PRIVATE
48!EM XIOS-OASIS-MCT compliance
49   PUBLIC   sbc_cpl_init       ! routine called by sbcmod.F90
50   PUBLIC   sbc_cpl_rcv        ! routine called by sbc_ice_lim(_2).F90
51   PUBLIC   sbc_cpl_snd        ! routine called by step.F90
52   PUBLIC   sbc_cpl_ice_tau    ! routine called by sbc_ice_lim(_2).F90
53   PUBLIC   sbc_cpl_ice_flx    ! routine called by sbc_ice_lim(_2).F90
54   PUBLIC   sbc_cpl_alloc      ! routine called in sbcice_cice.F90
55
56   INTEGER, PARAMETER ::   jpr_otx1   =  1            ! 3 atmosphere-ocean stress components on grid 1
57   INTEGER, PARAMETER ::   jpr_oty1   =  2            !
58   INTEGER, PARAMETER ::   jpr_otz1   =  3            !
59   INTEGER, PARAMETER ::   jpr_otx2   =  4            ! 3 atmosphere-ocean stress components on grid 2
60   INTEGER, PARAMETER ::   jpr_oty2   =  5            !
61   INTEGER, PARAMETER ::   jpr_otz2   =  6            !
62   INTEGER, PARAMETER ::   jpr_itx1   =  7            ! 3 atmosphere-ice   stress components on grid 1
63   INTEGER, PARAMETER ::   jpr_ity1   =  8            !
64   INTEGER, PARAMETER ::   jpr_itz1   =  9            !
65   INTEGER, PARAMETER ::   jpr_itx2   = 10            ! 3 atmosphere-ice   stress components on grid 2
66   INTEGER, PARAMETER ::   jpr_ity2   = 11            !
67   INTEGER, PARAMETER ::   jpr_itz2   = 12            !
68   INTEGER, PARAMETER ::   jpr_qsroce = 13            ! Qsr above the ocean
69   INTEGER, PARAMETER ::   jpr_qsrice = 14            ! Qsr above the ice
70   INTEGER, PARAMETER ::   jpr_qsrmix = 15 
71   INTEGER, PARAMETER ::   jpr_qnsoce = 16            ! Qns above the ocean
72   INTEGER, PARAMETER ::   jpr_qnsice = 17            ! Qns above the ice
73   INTEGER, PARAMETER ::   jpr_qnsmix = 18
74   INTEGER, PARAMETER ::   jpr_rain   = 19            ! total liquid precipitation (rain)
75   INTEGER, PARAMETER ::   jpr_snow   = 20            ! solid precipitation over the ocean (snow)
76   INTEGER, PARAMETER ::   jpr_tevp   = 21            ! total evaporation
77   INTEGER, PARAMETER ::   jpr_ievp   = 22            ! solid evaporation (sublimation)
78   INTEGER, PARAMETER ::   jpr_sbpr   = 23            ! sublimation - liquid precipitation - solid precipitation
79   INTEGER, PARAMETER ::   jpr_semp   = 24            ! solid freshwater budget (sublimation - snow)
80   INTEGER, PARAMETER ::   jpr_oemp   = 25            ! ocean freshwater budget (evap - precip)
81   INTEGER, PARAMETER ::   jpr_w10m   = 26            ! 10m wind
82   INTEGER, PARAMETER ::   jpr_dqnsdt = 27            ! d(Q non solar)/d(temperature)
83   INTEGER, PARAMETER ::   jpr_rnf    = 28            ! runoffs
84   INTEGER, PARAMETER ::   jpr_cal    = 29            ! calving
85   INTEGER, PARAMETER ::   jpr_taum   = 30            ! wind stress module
86   INTEGER, PARAMETER ::   jpr_co2    = 31
87   INTEGER, PARAMETER ::   jpr_topm   = 32            ! topmeltn
88   INTEGER, PARAMETER ::   jpr_botm   = 33            ! botmeltn
89   INTEGER, PARAMETER ::   jpr_ts_ice = 34            ! skin temperature of sea-ice (used for melt-ponds)
90   INTEGER, PARAMETER ::   jprcv      = 34            ! total number of fields received
91
92   INTEGER, PARAMETER ::   jps_fice   =  1            ! ice fraction
93   INTEGER, PARAMETER ::   jps_toce   =  2            ! ocean temperature
94   INTEGER, PARAMETER ::   jps_tice   =  3            ! ice   temperature
95   INTEGER, PARAMETER ::   jps_tmix   =  4            ! mixed temperature (ocean+ice)
96   INTEGER, PARAMETER ::   jps_albice =  5            ! ice   albedo
97   INTEGER, PARAMETER ::   jps_albmix =  6            ! mixed albedo
98   INTEGER, PARAMETER ::   jps_hice   =  7            ! ice  thickness
99   INTEGER, PARAMETER ::   jps_hsnw   =  8            ! snow thickness
100   INTEGER, PARAMETER ::   jps_ocx1   =  9            ! ocean current on grid 1
101   INTEGER, PARAMETER ::   jps_ocy1   = 10            !
102   INTEGER, PARAMETER ::   jps_ocz1   = 11            !
103   INTEGER, PARAMETER ::   jps_ivx1   = 12            ! ice   current on grid 1
104   INTEGER, PARAMETER ::   jps_ivy1   = 13            !
105   INTEGER, PARAMETER ::   jps_ivz1   = 14            !
106   INTEGER, PARAMETER ::   jps_co2    = 15
107   INTEGER, PARAMETER ::   jps_a_p    = 16            ! meltpond fraction 
108   INTEGER, PARAMETER ::   jps_ht_p   = 17            ! meltpond depth (m)
109   INTEGER, PARAMETER ::   jpsnd      = 18            ! total number of fields sent
110   !                                                         !!** namelist namsbc_cpl **
111   TYPE ::   FLD_C
112      CHARACTER(len = 32) ::   cldes                  ! desciption of the coupling strategy
113      CHARACTER(len = 32) ::   clcat                  ! multiple ice categories strategy
114      CHARACTER(len = 32) ::   clvref                 ! reference of vector ('spherical' or 'cartesian')
115      CHARACTER(len = 32) ::   clvor                  ! orientation of vector fields ('eastward-northward' or 'local grid')
116      CHARACTER(len = 32) ::   clvgrd                 ! grids on which is located the vector fields
117   END TYPE FLD_C
118   ! Send to the atmosphere                           !
119   TYPE(FLD_C) ::   sn_snd_temp, sn_snd_alb, sn_snd_thick, sn_snd_crt, sn_snd_co2, sn_snd_mpnd
120   ! Received from the atmosphere                     !
121   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
122   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_ts_ice 
123   ! Other namelist parameters                        !
124   INTEGER     ::   nn_cplmodel            ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
125   LOGICAL     ::   ln_usecplmask          !  use a coupling mask file to merge data received from several models
126                                           !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
127
128   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: xcplmask
129
130   TYPE ::   DYNARR     
131      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z3   
132   END TYPE DYNARR
133
134   TYPE( DYNARR ), SAVE, DIMENSION(jprcv) ::   frcv                      ! all fields recieved from the atmosphere
135
136   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   albedo_oce_mix     ! ocean albedo sent to atmosphere (mix clear/overcast sky)
137
138   INTEGER , ALLOCATABLE, SAVE, DIMENSION(    :) ::   nrcvinfo           ! OASIS info argument
139
140   !! Substitution
141#  include "vectopt_loop_substitute.h90"
142   !!----------------------------------------------------------------------
143   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
144   !! $Id$
145   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
146   !!----------------------------------------------------------------------
147
148CONTAINS
149 
150   INTEGER FUNCTION sbc_cpl_alloc()
151      !!----------------------------------------------------------------------
152      !!             ***  FUNCTION sbc_cpl_alloc  ***
153      !!----------------------------------------------------------------------
154      INTEGER :: ierr(3)
155      !!----------------------------------------------------------------------
156      ierr(:) = 0
157      !
158      ALLOCATE( albedo_oce_mix(jpi,jpj), nrcvinfo(jprcv),  STAT=ierr(1) )
159     
160#if ! defined key_lim3 && ! defined key_lim2 && ! defined key_cice
161      ALLOCATE( a_i(jpi,jpj,1) , STAT=ierr(2) )  ! used in sbcice_if.F90 (done here as there is no sbc_ice_if_init)
162#endif
163      ALLOCATE( xcplmask(jpi,jpj,nn_cplmodel) , STAT=ierr(3) )
164      !
165      sbc_cpl_alloc = MAXVAL( ierr )
166      IF( lk_mpp            )   CALL mpp_sum ( sbc_cpl_alloc )
167      IF( sbc_cpl_alloc > 0 )   CALL ctl_warn('sbc_cpl_alloc: allocation of arrays failed')
168      !
169   END FUNCTION sbc_cpl_alloc
170
171
172   SUBROUTINE sbc_cpl_init( k_ice )     
173      !!----------------------------------------------------------------------
174      !!             ***  ROUTINE sbc_cpl_init  ***
175      !!
176      !! ** Purpose :   Initialisation of send and received information from
177      !!                the atmospheric component
178      !!
179      !! ** Method  : * Read namsbc_cpl namelist
180      !!              * define the receive interface
181      !!              * define the send    interface
182      !!              * initialise the OASIS coupler
183      !!----------------------------------------------------------------------
184      INTEGER, INTENT(in) ::   k_ice    ! ice management in the sbc (=0/1/2/3)
185      !!
186      INTEGER ::   jn   ! dummy loop index
187      INTEGER ::   ios  ! Local integer output status for namelist read
188      INTEGER ::   inum 
189      REAL(wp), POINTER, DIMENSION(:,:) ::   zacs, zaos
190      !!
191      NAMELIST/namsbc_cpl/  sn_snd_temp, sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2,      &
192         &                  sn_snd_mpnd,                                                              &
193         &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr,      &
194         &                  sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf  , sn_rcv_cal   , sn_rcv_iceflx,   &
195         &                  sn_rcv_co2 , sn_rcv_ts_ice, nn_cplmodel , ln_usecplmask
196      !!---------------------------------------------------------------------
197      !
198      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_init')
199      !
200      CALL wrk_alloc( jpi,jpj, zacs, zaos )
201
202      ! ================================ !
203      !      Namelist informations       !
204      ! ================================ !
205
206      REWIND( numnam_ref )              ! Namelist namsbc_cpl in reference namelist : Variables for OASIS coupling
207      READ  ( numnam_ref, namsbc_cpl, IOSTAT = ios, ERR = 901)
208901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cpl in reference namelist', lwp )
209
210      REWIND( numnam_cfg )              ! Namelist namsbc_cpl in configuration namelist : Variables for OASIS coupling
211      READ  ( numnam_cfg, namsbc_cpl, IOSTAT = ios, ERR = 902 )
212902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cpl in configuration namelist', lwp )
213      IF(lwm) WRITE ( numond, namsbc_cpl )
214
215      IF(lwp) THEN                        ! control print
216         WRITE(numout,*)
217         WRITE(numout,*)'sbc_cpl_init : namsbc_cpl namelist '
218         WRITE(numout,*)'~~~~~~~~~~~~'
219         WRITE(numout,*)'  received fields (mutiple ice categogies)'
220         WRITE(numout,*)'      10m wind module                 = ', TRIM(sn_rcv_w10m%cldes  ), ' (', TRIM(sn_rcv_w10m%clcat  ), ')'
221         WRITE(numout,*)'      stress module                   = ', TRIM(sn_rcv_taumod%cldes), ' (', TRIM(sn_rcv_taumod%clcat), ')'
222         WRITE(numout,*)'      surface stress                  = ', TRIM(sn_rcv_tau%cldes   ), ' (', TRIM(sn_rcv_tau%clcat   ), ')'
223         WRITE(numout,*)'                     - referential    = ', sn_rcv_tau%clvref
224         WRITE(numout,*)'                     - orientation    = ', sn_rcv_tau%clvor
225         WRITE(numout,*)'                     - mesh           = ', sn_rcv_tau%clvgrd
226         WRITE(numout,*)'      non-solar heat flux sensitivity = ', TRIM(sn_rcv_dqnsdt%cldes), ' (', TRIM(sn_rcv_dqnsdt%clcat), ')'
227         WRITE(numout,*)'      solar heat flux                 = ', TRIM(sn_rcv_qsr%cldes   ), ' (', TRIM(sn_rcv_qsr%clcat   ), ')'
228         WRITE(numout,*)'      non-solar heat flux             = ', TRIM(sn_rcv_qns%cldes   ), ' (', TRIM(sn_rcv_qns%clcat   ), ')'
229         WRITE(numout,*)'      freshwater budget               = ', TRIM(sn_rcv_emp%cldes   ), ' (', TRIM(sn_rcv_emp%clcat   ), ')'
230         WRITE(numout,*)'      runoffs                         = ', TRIM(sn_rcv_rnf%cldes   ), ' (', TRIM(sn_rcv_rnf%clcat   ), ')'
231         WRITE(numout,*)'      calving                         = ', TRIM(sn_rcv_cal%cldes   ), ' (', TRIM(sn_rcv_cal%clcat   ), ')'
232         WRITE(numout,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')'
233         WRITE(numout,*)'      atm co2                         = ', TRIM(sn_rcv_co2%cldes   ), ' (', TRIM(sn_rcv_co2%clcat   ), ')'
234         WRITE(numout,*)'  sent fields (multiple ice categories)'
235         WRITE(numout,*)'      surface temperature             = ', TRIM(sn_snd_temp%cldes  ), ' (', TRIM(sn_snd_temp%clcat  ), ')'
236         WRITE(numout,*)'      albedo                          = ', TRIM(sn_snd_alb%cldes   ), ' (', TRIM(sn_snd_alb%clcat   ), ')'
237         WRITE(numout,*)'      ice/snow thickness              = ', TRIM(sn_snd_thick%cldes ), ' (', TRIM(sn_snd_thick%clcat ), ')'
238         WRITE(numout,*)'      surface current                 = ', TRIM(sn_snd_crt%cldes   ), ' (', TRIM(sn_snd_crt%clcat   ), ')'
239         WRITE(numout,*)'                      - referential   = ', sn_snd_crt%clvref 
240         WRITE(numout,*)'                      - orientation   = ', sn_snd_crt%clvor
241         WRITE(numout,*)'                      - mesh          = ', sn_snd_crt%clvgrd
242         WRITE(numout,*)'      oce co2 flux                    = ', TRIM(sn_snd_co2%cldes   ), ' (', TRIM(sn_snd_co2%clcat   ), ')'
243         WRITE(numout,*)'      meltponds fraction & depth      = ', TRIM(sn_snd_mpnd%cldes  ), ' (', TRIM(sn_snd_mpnd%clcat   ), ')'
244         WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel
245         WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask
246      ENDIF
247
248      !                                   ! allocate sbccpl arrays
249      IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )
250     
251      ! ================================ !
252      !   Define the receive interface   !
253      ! ================================ !
254      nrcvinfo(:) = OASIS_idle   ! needed by nrcvinfo(jpr_otx1) if we do not receive ocean stress
255
256      ! for each field: define the OASIS name                              (srcv(:)%clname)
257      !                 define receive or not from the namelist parameters (srcv(:)%laction)
258      !                 define the north fold type of lbc                  (srcv(:)%nsgn)
259
260      ! default definitions of srcv
261      srcv(:)%laction = .FALSE.   ;   srcv(:)%clgrid = 'T'   ;   srcv(:)%nsgn = 1.   ;   srcv(:)%nct = 1
262
263      !                                                      ! ------------------------- !
264      !                                                      ! ice and ocean wind stress !   
265      !                                                      ! ------------------------- !
266      !                                                           ! Name
267      srcv(jpr_otx1)%clname = 'O_OTaux1'      ! 1st ocean component on grid ONE (T or U)
268      srcv(jpr_oty1)%clname = 'O_OTauy1'      ! 2nd   -      -         -     -
269      srcv(jpr_otz1)%clname = 'O_OTauz1'      ! 3rd   -      -         -     -
270      srcv(jpr_otx2)%clname = 'O_OTaux2'      ! 1st ocean component on grid TWO (V)
271      srcv(jpr_oty2)%clname = 'O_OTauy2'      ! 2nd   -      -         -     -
272      srcv(jpr_otz2)%clname = 'O_OTauz2'      ! 3rd   -      -         -     -
273      !
274      srcv(jpr_itx1)%clname = 'O_ITaux1'      ! 1st  ice  component on grid ONE (T, F, I or U)
275      srcv(jpr_ity1)%clname = 'O_ITauy1'      ! 2nd   -      -         -     -
276      srcv(jpr_itz1)%clname = 'O_ITauz1'      ! 3rd   -      -         -     -
277      srcv(jpr_itx2)%clname = 'O_ITaux2'      ! 1st  ice  component on grid TWO (V)
278      srcv(jpr_ity2)%clname = 'O_ITauy2'      ! 2nd   -      -         -     -
279      srcv(jpr_itz2)%clname = 'O_ITauz2'      ! 3rd   -      -         -     -
280      !
281      ! Vectors: change of sign at north fold ONLY if on the local grid
282      IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' )   srcv(jpr_otx1:jpr_itz2)%nsgn = -1.
283     
284      !                                                           ! Set grid and action
285      SELECT CASE( TRIM( sn_rcv_tau%clvgrd ) )      !  'T', 'U,V', 'U,V,I', 'U,V,F', 'T,I', 'T,F', or 'T,U,V'
286      CASE( 'T' ) 
287         srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
288         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
289         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
290      CASE( 'U,V' ) 
291         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
292         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
293         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'U'        ! ice components given at U-point
294         srcv(jpr_itx2:jpr_itz2)%clgrid  = 'V'        !           and           V-point
295         srcv(jpr_otx1:jpr_itz2)%laction = .TRUE.     ! receive oce and ice components on both grid 1 & 2
296      CASE( 'U,V,T' )
297         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
298         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
299         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'T'        ! ice components given at T-point
300         srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
301         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
302      CASE( 'U,V,I' )
303         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
304         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
305         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'I'        ! ice components given at I-point
306         srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
307         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
308      CASE( 'U,V,F' )
309         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
310         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
311         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'F'        ! ice components given at F-point
312         srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
313         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
314      CASE( 'T,I' ) 
315         srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
316         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'I'        ! ice components given at I-point
317         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
318         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
319      CASE( 'T,F' ) 
320         srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
321         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'F'        ! ice components given at F-point
322         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
323         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
324      CASE( 'T,U,V' )
325         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'T'        ! oce components given at T-point
326         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'U'        ! ice components given at U-point
327         srcv(jpr_itx2:jpr_itz2)%clgrid  = 'V'        !           and           V-point
328         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1 only
329         srcv(jpr_itx1:jpr_itz2)%laction = .TRUE.     ! receive ice components on grid 1 & 2
330      CASE default   
331         CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_tau%clvgrd' )
332      END SELECT
333      !
334      IF( TRIM( sn_rcv_tau%clvref ) == 'spherical' )   &           ! spherical: 3rd component not received
335         &     srcv( (/jpr_otz1, jpr_otz2, jpr_itz1, jpr_itz2/) )%laction = .FALSE. 
336      !
337      IF( TRIM( sn_rcv_tau%clvor  ) == 'local grid' ) THEN        ! already on local grid -> no need of the second grid
338            srcv(jpr_otx2:jpr_otz2)%laction = .FALSE. 
339            srcv(jpr_itx2:jpr_itz2)%laction = .FALSE. 
340            srcv(jpr_oty1)%clgrid = srcv(jpr_oty2)%clgrid   ! not needed but cleaner...
341            srcv(jpr_ity1)%clgrid = srcv(jpr_ity2)%clgrid   ! not needed but cleaner...
342      ENDIF
343      !
344      IF( TRIM( sn_rcv_tau%cldes ) /= 'oce and ice' ) THEN        ! 'oce and ice' case ocean stress on ocean mesh used
345         srcv(jpr_itx1:jpr_itz2)%laction = .FALSE.    ! ice components not received
346         srcv(jpr_itx1)%clgrid = 'U'                  ! ocean stress used after its transformation
347         srcv(jpr_ity1)%clgrid = 'V'                  ! i.e. it is always at U- & V-points for i- & j-comp. resp.
348      ENDIF
349       
350      !                                                      ! ------------------------- !
351      !                                                      !    freshwater budget      !   E-P
352      !                                                      ! ------------------------- !
353      ! we suppose that atmosphere modele do not make the difference between precipiration (liquide or solid)
354      ! over ice of free ocean within the same atmospheric cell.cd
355      srcv(jpr_rain)%clname = 'OTotRain'      ! Rain = liquid precipitation
356      srcv(jpr_snow)%clname = 'OTotSnow'      ! Snow = solid precipitation
357      srcv(jpr_tevp)%clname = 'OTotEvap'      ! total evaporation (over oce + ice sublimation)
358      srcv(jpr_ievp)%clname = 'OIceEvap'      ! evaporation over ice = sublimation
359      srcv(jpr_sbpr)%clname = 'OSubMPre'      ! sublimation - liquid precipitation - solid precipitation
360      srcv(jpr_semp)%clname = 'OISubMSn'      ! ice solid water budget = sublimation - solid precipitation
361      srcv(jpr_oemp)%clname = 'OOEvaMPr'      ! ocean water budget = ocean Evap - ocean precip
362      SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
363      CASE( 'oce only'      )   ;   srcv(                                 jpr_oemp   )%laction = .TRUE. 
364      CASE( 'conservative'  )
365         srcv( (/jpr_rain, jpr_snow, jpr_ievp, jpr_tevp/) )%laction = .TRUE.
366         IF ( k_ice <= 1 )  srcv(jpr_ievp)%laction = .FALSE.
367      CASE( 'oce and ice'   )   ;   srcv( (/jpr_ievp, jpr_sbpr, jpr_semp, jpr_oemp/) )%laction = .TRUE.
368      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_emp%cldes' )
369      END SELECT
370      !Set the number of categories for coupling of sublimation
371      IF ( TRIM( sn_rcv_emp%clcat ) == 'yes' ) srcv(jpr_ievp)%nct = jpl
372      !
373      !                                                      ! ------------------------- !
374      !                                                      !     Runoffs & Calving     !   
375      !                                                      ! ------------------------- !
376      srcv(jpr_rnf   )%clname = 'O_Runoff'   ;   IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' )   srcv(jpr_rnf)%laction = .TRUE.
377! This isn't right - really just want ln_rnf_emp changed
378!                                                 IF( TRIM( sn_rcv_rnf%cldes ) == 'climato' )   THEN   ;   ln_rnf = .TRUE.
379!                                                 ELSE                                                 ;   ln_rnf = .FALSE.
380!                                                 ENDIF
381      srcv(jpr_cal   )%clname = 'OCalving'   ;   IF( TRIM( sn_rcv_cal%cldes ) == 'coupled' )   srcv(jpr_cal)%laction = .TRUE.
382
383      !                                                      ! ------------------------- !
384      !                                                      !    non solar radiation    !   Qns
385      !                                                      ! ------------------------- !
386      srcv(jpr_qnsoce)%clname = 'O_QnsOce'
387      srcv(jpr_qnsice)%clname = 'O_QnsIce'
388      srcv(jpr_qnsmix)%clname = 'O_QnsMix'
389      SELECT CASE( TRIM( sn_rcv_qns%cldes ) )
390      CASE( 'oce only'      )   ;   srcv(               jpr_qnsoce   )%laction = .TRUE.
391      CASE( 'conservative'  )   ;   srcv( (/jpr_qnsice, jpr_qnsmix/) )%laction = .TRUE.
392      CASE( 'oce and ice'   )   ;   srcv( (/jpr_qnsice, jpr_qnsoce/) )%laction = .TRUE.
393      CASE( 'mixed oce-ice' )   ;   srcv(               jpr_qnsmix   )%laction = .TRUE. 
394      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_qns%cldes' )
395      END SELECT
396      IF( TRIM( sn_rcv_qns%cldes ) == 'mixed oce-ice' .AND. jpl > 1 ) &
397         CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qns%cldes not currently allowed to be mixed oce-ice for multi-category ice' )
398      !                                                      ! ------------------------- !
399      !                                                      !    solar radiation        !   Qsr
400      !                                                      ! ------------------------- !
401      srcv(jpr_qsroce)%clname = 'O_QsrOce'
402      srcv(jpr_qsrice)%clname = 'O_QsrIce'
403      srcv(jpr_qsrmix)%clname = 'O_QsrMix'
404      SELECT CASE( TRIM( sn_rcv_qsr%cldes ) )
405      CASE( 'oce only'      )   ;   srcv(               jpr_qsroce   )%laction = .TRUE.
406      CASE( 'conservative'  )   ;   srcv( (/jpr_qsrice, jpr_qsrmix/) )%laction = .TRUE.
407      CASE( 'oce and ice'   )   ;   srcv( (/jpr_qsrice, jpr_qsroce/) )%laction = .TRUE.
408      CASE( 'mixed oce-ice' )   ;   srcv(               jpr_qsrmix   )%laction = .TRUE. 
409      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_qsr%cldes' )
410      END SELECT
411      IF( TRIM( sn_rcv_qsr%cldes ) == 'mixed oce-ice' .AND. jpl > 1 ) &
412         CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qsr%cldes not currently allowed to be mixed oce-ice for multi-category ice' )
413      !                                                      ! ------------------------- !
414      !                                                      !   non solar sensitivity   !   d(Qns)/d(T)
415      !                                                      ! ------------------------- !
416      srcv(jpr_dqnsdt)%clname = 'O_dQnsdT'   
417      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'coupled' )   srcv(jpr_dqnsdt)%laction = .TRUE.
418      !
419      ! non solar sensitivity mandatory for LIM ice model
420      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'none' .AND. k_ice /= 0 .AND. k_ice /= 4) &
421         CALL ctl_stop( 'sbc_cpl_init: sn_rcv_dqnsdt%cldes must be coupled in namsbc_cpl namelist' )
422      ! non solar sensitivity mandatory for mixed oce-ice solar radiation coupling technique
423      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'none' .AND. TRIM( sn_rcv_qns%cldes ) == 'mixed oce-ice' ) &
424         CALL ctl_stop( 'sbc_cpl_init: namsbc_cpl namelist mismatch between sn_rcv_qns%cldes and sn_rcv_dqnsdt%cldes' )
425      !                                                      ! ------------------------- !
426      !                                                      !      10m wind module      !   
427      !                                                      ! ------------------------- !
428      srcv(jpr_w10m)%clname = 'O_Wind10'   ;   IF( TRIM(sn_rcv_w10m%cldes  ) == 'coupled' )   srcv(jpr_w10m)%laction = .TRUE. 
429      !
430      !                                                      ! ------------------------- !
431      !                                                      !   wind stress module      !   
432      !                                                      ! ------------------------- !
433      srcv(jpr_taum)%clname = 'O_TauMod'   ;   IF( TRIM(sn_rcv_taumod%cldes) == 'coupled' )   srcv(jpr_taum)%laction = .TRUE.
434      lhftau = srcv(jpr_taum)%laction
435
436      !                                                      ! ------------------------- !
437      !                                                      !      Atmospheric CO2      !
438      !                                                      ! ------------------------- !
439      srcv(jpr_co2 )%clname = 'O_AtmCO2'   ;   IF( TRIM(sn_rcv_co2%cldes   ) == 'coupled' )    srcv(jpr_co2 )%laction = .TRUE.
440      !                                                      ! ------------------------- !
441      !                                                      !   topmelt and botmelt     !   
442      !                                                      ! ------------------------- !
443      srcv(jpr_topm )%clname = 'OTopMlt'
444      srcv(jpr_botm )%clname = 'OBotMlt'
445      IF( TRIM(sn_rcv_iceflx%cldes) == 'coupled' ) THEN
446         IF ( TRIM( sn_rcv_iceflx%clcat ) == 'yes' ) THEN
447            srcv(jpr_topm:jpr_botm)%nct = jpl
448         ELSE
449            CALL ctl_stop( 'sbc_cpl_init: sn_rcv_iceflx%clcat should always be set to yes currently' )
450         ENDIF
451         srcv(jpr_topm:jpr_botm)%laction = .TRUE.
452      ENDIF
453     
454#if defined key_cice && ! defined key_cice4
455      !                                                      ! ----------------------------- !
456      !                                                      !  sea-ice skin temperature     !   
457      !                                                      !  used in meltpond scheme      !
458      !                                                      !  May be calculated in Atm     !
459      !                                                      ! ----------------------------- !
460      srcv(jpr_ts_ice)%clname = 'OTsfIce'
461      IF ( TRIM( sn_rcv_ts_ice%cldes ) == 'ice' ) srcv(jpr_ts_ice)%laction = .TRUE.
462      IF ( TRIM( sn_rcv_ts_ice%clcat ) == 'yes' ) srcv(jpr_ts_ice)%nct = jpl
463      !TODO: Should there be a consistency check here?
464#endif
465
466      ! Allocate all parts of frcv used for received fields
467      DO jn = 1, jprcv
468         IF ( srcv(jn)%laction ) ALLOCATE( frcv(jn)%z3(jpi,jpj,srcv(jn)%nct) )
469      END DO
470      ! Allocate taum part of frcv which is used even when not received as coupling field
471      IF ( .NOT. srcv(jpr_taum)%laction ) ALLOCATE( frcv(jpr_taum)%z3(jpi,jpj,srcv(jpr_taum)%nct) )
472      ! Allocate itx1 and ity1 as they are used in sbc_cpl_ice_tau even if srcv(jpr_itx1)%laction = .FALSE.
473      IF( k_ice /= 0 ) THEN
474         IF ( .NOT. srcv(jpr_itx1)%laction ) ALLOCATE( frcv(jpr_itx1)%z3(jpi,jpj,srcv(jpr_itx1)%nct) )
475         IF ( .NOT. srcv(jpr_ity1)%laction ) ALLOCATE( frcv(jpr_ity1)%z3(jpi,jpj,srcv(jpr_ity1)%nct) )
476      END IF
477
478      ! ================================ !
479      !     Define the send interface    !
480      ! ================================ !
481      ! for each field: define the OASIS name                           (ssnd(:)%clname)
482      !                 define send or not from the namelist parameters (ssnd(:)%laction)
483      !                 define the north fold type of lbc               (ssnd(:)%nsgn)
484     
485      ! default definitions of nsnd
486      ssnd(:)%laction = .FALSE.   ;   ssnd(:)%clgrid = 'T'   ;   ssnd(:)%nsgn = 1.  ; ssnd(:)%nct = 1
487         
488      !                                                      ! ------------------------- !
489      !                                                      !    Surface temperature    !
490      !                                                      ! ------------------------- !
491      ssnd(jps_toce)%clname = 'O_SSTSST'
492      ssnd(jps_tice)%clname = 'O_TepIce'
493      ssnd(jps_tmix)%clname = 'O_TepMix'
494      SELECT CASE( TRIM( sn_snd_temp%cldes ) )
495      CASE( 'none'         )       ! nothing to do
496      CASE( 'oce only'             )   ;   ssnd(   jps_toce             )%laction = .TRUE.
497      CASE( 'weighted oce and ice' )
498         ssnd( (/jps_toce, jps_tice/) )%laction = .TRUE.
499         IF ( TRIM( sn_snd_temp%clcat ) == 'yes' )  ssnd(jps_tice)%nct = jpl
500      CASE( 'mixed oce-ice'        )   ;   ssnd(   jps_tmix             )%laction = .TRUE.
501      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_temp%cldes' )
502      END SELECT
503     
504      !                                                      ! ------------------------- !
505      !                                                      !          Albedo           !
506      !                                                      ! ------------------------- !
507      ssnd(jps_albice)%clname = 'O_AlbIce' 
508      ssnd(jps_albmix)%clname = 'O_AlbMix'
509      SELECT CASE( TRIM( sn_snd_alb%cldes ) )
510      CASE( 'none'          )       ! nothing to do
511      CASE( 'weighted ice'  )   ;   ssnd(jps_albice)%laction = .TRUE.
512      CASE( 'mixed oce-ice' )   ;   ssnd(jps_albmix)%laction = .TRUE.
513      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_alb%cldes' )
514      END SELECT
515      !
516      ! Need to calculate oceanic albedo if
517      !     1. sending mixed oce-ice albedo or
518      !     2. receiving mixed oce-ice solar radiation
519      IF ( TRIM ( sn_snd_alb%cldes ) == 'mixed oce-ice' .OR. TRIM ( sn_rcv_qsr%cldes ) == 'mixed oce-ice' ) THEN
520         CALL albedo_oce( zaos, zacs )
521         ! Due to lack of information on nebulosity : mean clear/overcast sky
522         albedo_oce_mix(:,:) = ( zacs(:,:) + zaos(:,:) ) * 0.5
523      ENDIF
524
525      !                                                      ! ------------------------- !
526      !                                                      !  Ice fraction & Thickness
527      !                                                      ! ------------------------- !
528      ssnd(jps_fice)%clname = 'OIceFrc'
529      ssnd(jps_hice)%clname = 'OIceTck'
530      ssnd(jps_hsnw)%clname = 'OSnwTck'
531      ssnd(jps_a_p)%clname  = 'OPndFrc'
532      ssnd(jps_ht_p)%clname = 'OPndTck'
533      IF( k_ice /= 0 ) THEN
534         ssnd(jps_fice)%laction = .TRUE.                  ! if ice treated in the ocean (even in climato case)
535! Currently no namelist entry to determine sending of multi-category ice fraction so use the thickness entry for now
536         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_fice)%nct = jpl
537      ENDIF
538
539      SELECT CASE ( TRIM( sn_snd_thick%cldes ) )
540      CASE( 'none'         )       ! nothing to do
541      CASE( 'ice and snow' ) 
542         ssnd(jps_hice:jps_hsnw)%laction = .TRUE.
543         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) THEN
544            ssnd(jps_hice:jps_hsnw)%nct = jpl
545         ELSE
546            IF ( jpl > 1 ) THEN
547CALL ctl_stop( 'sbc_cpl_init: use weighted ice and snow option for sn_snd_thick%cldes if not exchanging category fields' )
548            ENDIF
549         ENDIF
550      CASE ( 'weighted ice and snow' ) 
551         ssnd(jps_hice:jps_hsnw)%laction = .TRUE.
552         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_hice:jps_hsnw)%nct = jpl
553      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_thick%cldes' )
554      END SELECT
555
556      !                                                      ! ------------------------- !
557      !                                                      ! Ice Meltponds             !
558      !                                                      ! ------------------------- !
559#if defined key_cice && ! defined key_cice4
560      ! Meltponds only CICE5
561      ssnd(jps_a_p)%clname = 'OPndFrc'   
562      ssnd(jps_ht_p)%clname = 'OPndTck'   
563      SELECT CASE ( TRIM( sn_snd_mpnd%cldes ) )
564      CASE ( 'none' )
565         ssnd(jps_a_p)%laction = .FALSE.
566         ssnd(jps_ht_p)%laction = .FALSE.
567      CASE ( 'ice only' ) 
568         ssnd(jps_a_p)%laction = .TRUE.
569         ssnd(jps_ht_p)%laction = .TRUE.
570         IF ( TRIM( sn_snd_mpnd%clcat ) == 'yes' ) THEN
571            ssnd(jps_a_p)%nct = jpl
572            ssnd(jps_ht_p)%nct = jpl
573         ELSE
574            IF ( jpl > 1 ) THEN
575               CALL ctl_stop( 'sbc_cpl_init: use weighted ice option for sn_snd_mpnd%cldes if not exchanging category fields' )
576            ENDIF
577         ENDIF
578      CASE ( 'weighted ice' ) 
579         ssnd(jps_a_p)%laction = .TRUE.
580         ssnd(jps_ht_p)%laction = .TRUE.
581         IF ( TRIM( sn_snd_mpnd%clcat ) == 'yes' ) THEN
582            ssnd(jps_a_p)%nct = jpl 
583            ssnd(jps_ht_p)%nct = jpl 
584         ENDIF
585      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_mpnd%cldes' )
586      END SELECT
587#else
588      IF( TRIM( sn_snd_mpnd%cldes /= 'none' ) THEN
589         CALL ctl_stop('Meltponds can only be used with CICEv5')
590      ENDIF
591#endif
592
593      !                                                      ! ------------------------- !
594      !                                                      !      Surface current      !
595      !                                                      ! ------------------------- !
596      !        ocean currents              !            ice velocities
597      ssnd(jps_ocx1)%clname = 'O_OCurx1'   ;   ssnd(jps_ivx1)%clname = 'O_IVelx1'
598      ssnd(jps_ocy1)%clname = 'O_OCury1'   ;   ssnd(jps_ivy1)%clname = 'O_IVely1'
599      ssnd(jps_ocz1)%clname = 'O_OCurz1'   ;   ssnd(jps_ivz1)%clname = 'O_IVelz1'
600      !
601      ssnd(jps_ocx1:jps_ivz1)%nsgn = -1.   ! vectors: change of the sign at the north fold
602
603      IF( sn_snd_crt%clvgrd == 'U,V' ) THEN
604         ssnd(jps_ocx1)%clgrid = 'U' ; ssnd(jps_ocy1)%clgrid = 'V'
605      ELSE IF( sn_snd_crt%clvgrd /= 'T' ) THEN 
606         CALL ctl_stop( 'sn_snd_crt%clvgrd must be equal to T' )
607         ssnd(jps_ocx1:jps_ivz1)%clgrid  = 'T'      ! all oce and ice components on the same unique grid
608      ENDIF
609      ssnd(jps_ocx1:jps_ivz1)%laction = .TRUE.   ! default: all are send
610      IF( TRIM( sn_snd_crt%clvref ) == 'spherical' )   ssnd( (/jps_ocz1, jps_ivz1/) )%laction = .FALSE. 
611      IF( TRIM( sn_snd_crt%clvor ) == 'eastward-northward' ) ssnd(jps_ocx1:jps_ivz1)%nsgn = 1.
612      SELECT CASE( TRIM( sn_snd_crt%cldes ) )
613      CASE( 'none'                 )   ;   ssnd(jps_ocx1:jps_ivz1)%laction = .FALSE.
614      CASE( 'oce only'             )   ;   ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE.
615      CASE( 'weighted oce and ice' )   !   nothing to do
616      CASE( 'mixed oce-ice'        )   ;   ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE.
617      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_crt%cldes' )
618      END SELECT
619
620      !                                                      ! ------------------------- !
621      !                                                      !          CO2 flux         !
622      !                                                      ! ------------------------- !
623      ssnd(jps_co2)%clname = 'O_CO2FLX' ;  IF( TRIM(sn_snd_co2%cldes) == 'coupled' )    ssnd(jps_co2 )%laction = .TRUE.
624      !
625      ! ================================ !
626      !   initialisation of the coupler  !
627      ! ================================ !
628
629      CALL cpl_define(jprcv, jpsnd,nn_cplmodel)           
630      IF (ln_usecplmask) THEN
631         xcplmask(:,:,:) = 0.
632         CALL iom_open( 'cplmask', inum )
633         CALL iom_get( inum, jpdom_unknown, 'cplmask', xcplmask(1:nlci,1:nlcj,1:nn_cplmodel),   &
634            &          kstart = (/ mig(1),mjg(1),1 /), kcount = (/ nlci,nlcj,nn_cplmodel /) )
635         CALL iom_close( inum )
636      ELSE
637         xcplmask(:,:,:) = 1.
638      ENDIF
639      !
640      IF( ln_dm2dc .AND. ( cpl_freq( jpr_qsroce ) + cpl_freq( jpr_qsrmix ) /= 86400 ) )   &
641         &   CALL ctl_stop( 'sbc_cpl_init: diurnal cycle reconstruction (ln_dm2dc) needs daily couping for solar radiation' )
642
643      CALL wrk_dealloc( jpi,jpj, zacs, zaos )
644      !
645      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_init')
646      !
647   END SUBROUTINE sbc_cpl_init
648
649
650   SUBROUTINE sbc_cpl_rcv( kt, k_fsbc, k_ice )     
651      !!----------------------------------------------------------------------
652      !!             ***  ROUTINE sbc_cpl_rcv  ***
653      !!
654      !! ** Purpose :   provide the stress over the ocean and, if no sea-ice,
655      !!                provide the ocean heat and freshwater fluxes.
656      !!
657      !! ** Method  : - Receive all the atmospheric fields (stored in frcv array). called at each time step.
658      !!                OASIS controls if there is something do receive or not. nrcvinfo contains the info
659      !!                to know if the field was really received or not
660      !!
661      !!              --> If ocean stress was really received:
662      !!
663      !!                  - transform the received ocean stress vector from the received
664      !!                 referential and grid into an atmosphere-ocean stress in
665      !!                 the (i,j) ocean referencial and at the ocean velocity point.
666      !!                    The received stress are :
667      !!                     - defined by 3 components (if cartesian coordinate)
668      !!                            or by 2 components (if spherical)
669      !!                     - oriented along geographical   coordinate (if eastward-northward)
670      !!                            or  along the local grid coordinate (if local grid)
671      !!                     - given at U- and V-point, resp.   if received on 2 grids
672      !!                            or at T-point               if received on 1 grid
673      !!                    Therefore and if necessary, they are successively
674      !!                  processed in order to obtain them
675      !!                     first  as  2 components on the sphere
676      !!                     second as  2 components oriented along the local grid
677      !!                     third  as  2 components on the U,V grid
678      !!
679      !!              -->
680      !!
681      !!              - In 'ocean only' case, non solar and solar ocean heat fluxes
682      !!             and total ocean freshwater fluxes 
683      !!
684      !! ** Method  :   receive all fields from the atmosphere and transform
685      !!              them into ocean surface boundary condition fields
686      !!
687      !! ** Action  :   update  utau, vtau   ocean stress at U,V grid
688      !!                        taum         wind stress module at T-point
689      !!                        wndm         wind speed  module at T-point over free ocean or leads in presence of sea-ice
690      !!                        qns          non solar heat fluxes including emp heat content    (ocean only case)
691      !!                                     and the latent heat flux of solid precip. melting
692      !!                        qsr          solar ocean heat fluxes   (ocean only case)
693      !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case)
694      !!----------------------------------------------------------------------
695      INTEGER, INTENT(in) ::   kt       ! ocean model time step index
696      INTEGER, INTENT(in) ::   k_fsbc   ! frequency of sbc (-> ice model) computation
697      INTEGER, INTENT(in) ::   k_ice    ! ice management in the sbc (=0/1/2/3)
698      !!
699      LOGICAL ::    llnewtx, llnewtau      ! update wind stress components and module??
700      INTEGER  ::   ji, jj, jl, jn         ! dummy loop indices
701      INTEGER  ::   isec                   ! number of seconds since nit000 (assuming rdttra did not change since nit000)
702      REAL(wp) ::   zcumulneg, zcumulpos   ! temporary scalars     
703      REAL(wp) ::   zcoef                  ! temporary scalar
704      REAL(wp) ::   zrhoa  = 1.22          ! Air density kg/m3
705      REAL(wp) ::   zcdrag = 1.5e-3        ! drag coefficient
706      REAL(wp) ::   zzx, zzy               ! temporary variables
707      REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty 
708      !!----------------------------------------------------------------------
709      !
710      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_rcv')
711      !
712      CALL wrk_alloc( jpi,jpj, ztx, zty )
713      !                                                 ! Receive all the atmos. fields (including ice information)
714      isec = ( kt - nit000 ) * NINT( rdttra(1) )             ! date of exchanges
715      DO jn = 1, jprcv                                       ! received fields sent by the atmosphere
716         IF( srcv(jn)%laction )   CALL cpl_rcv( jn, isec, frcv(jn)%z3, xcplmask, nrcvinfo(jn) )
717      END DO
718
719      !                                                      ! ========================= !
720      IF( srcv(jpr_otx1)%laction ) THEN                      !  ocean stress components  !
721         !                                                   ! ========================= !
722         ! define frcv(jpr_otx1)%z3(:,:,1) and frcv(jpr_oty1)%z3(:,:,1): stress at U/V point along model grid
723         ! => need to be done only when we receive the field
724         IF(  nrcvinfo(jpr_otx1) == OASIS_Rcv ) THEN
725            !
726            IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN            ! 2 components on the sphere
727               !                                                       ! (cartesian to spherical -> 3 to 2 components)
728               !
729               CALL geo2oce( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), frcv(jpr_otz1)%z3(:,:,1),   &
730                  &          srcv(jpr_otx1)%clgrid, ztx, zty )
731               frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 1st grid
732               frcv(jpr_oty1)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 1st grid
733               !
734               IF( srcv(jpr_otx2)%laction ) THEN
735                  CALL geo2oce( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), frcv(jpr_otz2)%z3(:,:,1),   &
736                     &          srcv(jpr_otx2)%clgrid, ztx, zty )
737                  frcv(jpr_otx2)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 2nd grid
738                  frcv(jpr_oty2)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 2nd grid
739               ENDIF
740               !
741            ENDIF
742            !
743            IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN   ! 2 components oriented along the local grid
744               !                                                       ! (geographical to local grid -> rotate the components)
745               CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->i', ztx )   
746               IF( srcv(jpr_otx2)%laction ) THEN
747                  CALL rot_rep( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), srcv(jpr_otx2)%clgrid, 'en->j', zty )   
748               ELSE 
749                  CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->j', zty ) 
750               ENDIF
751               frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid
752               frcv(jpr_oty1)%z3(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 2nd grid
753            ENDIF
754            !                             
755            IF( srcv(jpr_otx1)%clgrid == 'T' ) THEN
756               DO jj = 2, jpjm1                                          ! T ==> (U,V)
757                  DO ji = fs_2, fs_jpim1   ! vector opt.
758                     frcv(jpr_otx1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_otx1)%z3(ji+1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1) )
759                     frcv(jpr_oty1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_oty1)%z3(ji  ,jj+1,1) + frcv(jpr_oty1)%z3(ji,jj,1) )
760                  END DO
761               END DO
762               CALL lbc_lnk( frcv(jpr_otx1)%z3(:,:,1), 'U',  -1. )   ;   CALL lbc_lnk( frcv(jpr_oty1)%z3(:,:,1), 'V',  -1. )
763            ENDIF
764            llnewtx = .TRUE.
765         ELSE
766            llnewtx = .FALSE.
767         ENDIF
768         !                                                   ! ========================= !
769      ELSE                                                   !   No dynamical coupling   !
770         !                                                   ! ========================= !
771         frcv(jpr_otx1)%z3(:,:,1) = 0.e0                               ! here simply set to zero
772         frcv(jpr_oty1)%z3(:,:,1) = 0.e0                               ! an external read in a file can be added instead
773         llnewtx = .TRUE.
774         !
775      ENDIF
776     
777      !                                                      ! ========================= !
778      !                                                      !    wind stress module     !   (taum)
779      !                                                      ! ========================= !
780      !
781      IF( .NOT. srcv(jpr_taum)%laction ) THEN                    ! compute wind stress module from its components if not received
782         ! => need to be done only when otx1 was changed
783         IF( llnewtx ) THEN
784!CDIR NOVERRCHK
785            DO jj = 2, jpjm1
786!CDIR NOVERRCHK
787               DO ji = fs_2, fs_jpim1   ! vect. opt.
788                  zzx = frcv(jpr_otx1)%z3(ji-1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1)
789                  zzy = frcv(jpr_oty1)%z3(ji  ,jj-1,1) + frcv(jpr_oty1)%z3(ji,jj,1)
790                  frcv(jpr_taum)%z3(ji,jj,1) = 0.5 * SQRT( zzx * zzx + zzy * zzy )
791               END DO
792            END DO
793            CALL lbc_lnk( frcv(jpr_taum)%z3(:,:,1), 'T', 1. )
794            llnewtau = .TRUE.
795         ELSE
796            llnewtau = .FALSE.
797         ENDIF
798      ELSE
799         llnewtau = nrcvinfo(jpr_taum) == OASIS_Rcv
800         ! Stress module can be negative when received (interpolation problem)
801         IF( llnewtau ) THEN
802            frcv(jpr_taum)%z3(:,:,1) = MAX( 0._wp, frcv(jpr_taum)%z3(:,:,1) )
803         ENDIF
804      ENDIF
805     
806      !                                                      ! ========================= !
807      !                                                      !      10 m wind speed      !   (wndm)
808      !                                                      ! ========================= !
809      !
810      IF( .NOT. srcv(jpr_w10m)%laction ) THEN                    ! compute wind spreed from wind stress module if not received 
811         ! => need to be done only when taumod was changed
812         IF( llnewtau ) THEN
813            zcoef = 1. / ( zrhoa * zcdrag ) 
814!CDIR NOVERRCHK
815            DO jj = 1, jpj
816!CDIR NOVERRCHK
817               DO ji = 1, jpi 
818                  wndm(ji,jj) = SQRT( frcv(jpr_taum)%z3(ji,jj,1) * zcoef )
819               END DO
820            END DO
821         ENDIF
822      ELSE
823         IF ( nrcvinfo(jpr_w10m) == OASIS_Rcv ) wndm(:,:) = frcv(jpr_w10m)%z3(:,:,1)
824      ENDIF
825
826      ! u(v)tau and taum will be modified by ice model
827      ! -> need to be reset before each call of the ice/fsbc     
828      IF( MOD( kt-1, k_fsbc ) == 0 ) THEN
829         !
830         utau(:,:) = frcv(jpr_otx1)%z3(:,:,1)
831         vtau(:,:) = frcv(jpr_oty1)%z3(:,:,1)
832         taum(:,:) = frcv(jpr_taum)%z3(:,:,1)
833         CALL iom_put( "taum_oce", taum )   ! output wind stress module
834         
835      ENDIF
836
837#if defined key_cpl_carbon_cycle
838      !                                                              ! atmosph. CO2 (ppm)
839      IF( srcv(jpr_co2)%laction )   atm_co2(:,:) = frcv(jpr_co2)%z3(:,:,1)
840#endif
841
842
843#if defined key_cice && ! defined key_cice4
844      !  ! Sea ice surface skin temp:
845      ! Use tn_ice for this as it's not used for anything else in cice case
846      IF( srcv(jpr_ts_ice)%laction ) THEN
847        DO jl = 1, jpl
848          DO jj = 1, jpj
849            DO ji = 1, jpi
850              IF (frcv(jpr_ts_ice)%z3(ji,jj,jl) > 0.0) THEN
851                tn_ice(ji,jj,jl) = 0.0
852              ELSE IF (frcv(jpr_ts_ice)%z3(ji,jj,jl) < -60.0) THEN
853                tn_ice(ji,jj,jl) = -60.0
854              ELSE
855                tn_ice(ji,jj,jl) = frcv(jpr_ts_ice)%z3(ji,jj,jl)
856              ENDIF
857            END DO
858          END DO
859        END DO
860      ENDIF
861#endif
862      !                                                      ! ========================= !
863      IF( k_ice <= 1 ) THEN                                  !  heat & freshwater fluxes ! (Ocean only case)
864         !                                                   ! ========================= !
865         !
866         !                                                       ! total freshwater fluxes over the ocean (emp)
867         SELECT CASE( TRIM( sn_rcv_emp%cldes ) )                                    ! evaporation - precipitation
868         CASE( 'conservative' )
869            emp(:,:) = frcv(jpr_tevp)%z3(:,:,1) - ( frcv(jpr_rain)%z3(:,:,1) + frcv(jpr_snow)%z3(:,:,1) )
870         CASE( 'oce only', 'oce and ice' )
871            emp(:,:) = frcv(jpr_oemp)%z3(:,:,1)
872         CASE default
873            CALL ctl_stop( 'sbc_cpl_rcv: wrong definition of sn_rcv_emp%cldes' )
874         END SELECT
875         !
876         !                                                        ! runoffs and calving (added in emp)
877         IF( srcv(jpr_rnf)%laction )   emp(:,:) = emp(:,:) - frcv(jpr_rnf)%z3(:,:,1)
878         IF( srcv(jpr_cal)%laction )   emp(:,:) = emp(:,:) - frcv(jpr_cal)%z3(:,:,1)
879         !
880!!gm :  this seems to be internal cooking, not sure to need that in a generic interface
881!!gm                                       at least should be optional...
882!!         IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' ) THEN     ! add to the total freshwater budget
883!!            ! remove negative runoff
884!!            zcumulpos = SUM( MAX( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * e1t(:,:) * e2t(:,:) * tmask_i(:,:) )
885!!            zcumulneg = SUM( MIN( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * e1t(:,:) * e2t(:,:) * tmask_i(:,:) )
886!!            IF( lk_mpp )   CALL mpp_sum( zcumulpos )   ! sum over the global domain
887!!            IF( lk_mpp )   CALL mpp_sum( zcumulneg )
888!!            IF( zcumulpos /= 0. ) THEN                 ! distribute negative runoff on positive runoff grid points
889!!               zcumulneg = 1.e0 + zcumulneg / zcumulpos
890!!               frcv(jpr_rnf)%z3(:,:,1) = MAX( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * zcumulneg
891!!            ENDIF     
892!!            ! add runoff to e-p
893!!            emp(:,:) = emp(:,:) - frcv(jpr_rnf)%z3(:,:,1)
894!!         ENDIF
895!!gm  end of internal cooking
896         !
897         !                                                       ! non solar heat flux over the ocean (qns)
898         IF( srcv(jpr_qnsoce)%laction )   qns(:,:) = frcv(jpr_qnsoce)%z3(:,:,1)
899         IF( srcv(jpr_qnsmix)%laction )   qns(:,:) = frcv(jpr_qnsmix)%z3(:,:,1)
900         ! update qns over the free ocean with:
901         qns(:,:) =  qns(:,:) - emp(:,:) * sst_m(:,:) * rcp            ! remove heat content due to mass flux (assumed to be at SST)
902         IF( srcv(jpr_snow  )%laction )   THEN
903              qns(:,:) = qns(:,:) - frcv(jpr_snow)%z3(:,:,1) * lfus    ! energy for melting solid precipitation over the free ocean
904         ENDIF
905
906         !                                                       ! solar flux over the ocean          (qsr)
907         IF( srcv(jpr_qsroce)%laction )   qsr(:,:) = frcv(jpr_qsroce)%z3(:,:,1)
908         IF( srcv(jpr_qsrmix)%laction )   qsr(:,:) = frcv(jpr_qsrmix)%z3(:,:,1)
909         IF( ln_dm2dc )   qsr(:,:) = sbc_dcy( qsr )                           ! modify qsr to include the diurnal cycle
910         !
911 
912      ENDIF
913      !
914      CALL wrk_dealloc( jpi,jpj, ztx, zty )
915      !
916      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_rcv')
917      !
918   END SUBROUTINE sbc_cpl_rcv
919   
920
921   SUBROUTINE sbc_cpl_ice_tau( p_taui, p_tauj )     
922      !!----------------------------------------------------------------------
923      !!             ***  ROUTINE sbc_cpl_ice_tau  ***
924      !!
925      !! ** Purpose :   provide the stress over sea-ice in coupled mode
926      !!
927      !! ** Method  :   transform the received stress from the atmosphere into
928      !!             an atmosphere-ice stress in the (i,j) ocean referencial
929      !!             and at the velocity point of the sea-ice model (cp_ice_msh):
930      !!                'C'-grid : i- (j-) components given at U- (V-) point
931      !!                'I'-grid : B-grid lower-left corner: both components given at I-point
932      !!
933      !!                The received stress are :
934      !!                 - defined by 3 components (if cartesian coordinate)
935      !!                        or by 2 components (if spherical)
936      !!                 - oriented along geographical   coordinate (if eastward-northward)
937      !!                        or  along the local grid coordinate (if local grid)
938      !!                 - given at U- and V-point, resp.   if received on 2 grids
939      !!                        or at a same point (T or I) if received on 1 grid
940      !!                Therefore and if necessary, they are successively
941      !!             processed in order to obtain them
942      !!                 first  as  2 components on the sphere
943      !!                 second as  2 components oriented along the local grid
944      !!                 third  as  2 components on the cp_ice_msh point
945      !!
946      !!                Except in 'oce and ice' case, only one vector stress field
947      !!             is received. It has already been processed in sbc_cpl_rcv
948      !!             so that it is now defined as (i,j) components given at U-
949      !!             and V-points, respectively. Therefore, only the third
950      !!             transformation is done and only if the ice-grid is a 'I'-grid.
951      !!
952      !! ** Action  :   return ptau_i, ptau_j, the stress over the ice at cp_ice_msh point
953      !!----------------------------------------------------------------------
954      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_taui   ! i- & j-components of atmos-ice stress [N/m2]
955      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_tauj   ! at I-point (B-grid) or U & V-point (C-grid)
956      !!
957      INTEGER ::   ji, jj                          ! dummy loop indices
958      INTEGER ::   itx                             ! index of taux over ice
959      REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty 
960      !!----------------------------------------------------------------------
961      !
962      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_tau')
963      !
964      CALL wrk_alloc( jpi,jpj, ztx, zty )
965
966      IF( srcv(jpr_itx1)%laction ) THEN   ;   itx =  jpr_itx1   
967      ELSE                                ;   itx =  jpr_otx1
968      ENDIF
969
970      ! do something only if we just received the stress from atmosphere
971      IF(  nrcvinfo(itx) == OASIS_Rcv ) THEN
972
973         !                                                      ! ======================= !
974         IF( srcv(jpr_itx1)%laction ) THEN                      !   ice stress received   !
975            !                                                   ! ======================= !
976           
977            IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN            ! 2 components on the sphere
978               !                                                       ! (cartesian to spherical -> 3 to 2 components)
979               CALL geo2oce(  frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), frcv(jpr_itz1)%z3(:,:,1),   &
980                  &          srcv(jpr_itx1)%clgrid, ztx, zty )
981               frcv(jpr_itx1)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 1st grid
982               frcv(jpr_ity1)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 1st grid
983               !
984               IF( srcv(jpr_itx2)%laction ) THEN
985                  CALL geo2oce( frcv(jpr_itx2)%z3(:,:,1), frcv(jpr_ity2)%z3(:,:,1), frcv(jpr_itz2)%z3(:,:,1),   &
986                     &          srcv(jpr_itx2)%clgrid, ztx, zty )
987                  frcv(jpr_itx2)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 2nd grid
988                  frcv(jpr_ity2)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 2nd grid
989               ENDIF
990               !
991            ENDIF
992            !
993            IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN   ! 2 components oriented along the local grid
994               !                                                       ! (geographical to local grid -> rotate the components)
995               CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), srcv(jpr_itx1)%clgrid, 'en->i', ztx )   
996               IF( srcv(jpr_itx2)%laction ) THEN
997                  CALL rot_rep( frcv(jpr_itx2)%z3(:,:,1), frcv(jpr_ity2)%z3(:,:,1), srcv(jpr_itx2)%clgrid, 'en->j', zty )   
998               ELSE
999                  CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), srcv(jpr_itx1)%clgrid, 'en->j', zty ) 
1000               ENDIF
1001               frcv(jpr_itx1)%z3(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid
1002               frcv(jpr_ity1)%z3(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 1st grid
1003            ENDIF
1004            !                                                   ! ======================= !
1005         ELSE                                                   !     use ocean stress    !
1006            !                                                   ! ======================= !
1007            frcv(jpr_itx1)%z3(:,:,1) = frcv(jpr_otx1)%z3(:,:,1)
1008            frcv(jpr_ity1)%z3(:,:,1) = frcv(jpr_oty1)%z3(:,:,1)
1009            !
1010         ENDIF
1011
1012         !                                                      ! ======================= !
1013         !                                                      !     put on ice grid     !
1014         !                                                      ! ======================= !
1015         !   
1016         !                                                  j+1   j     -----V---F
1017         ! ice stress on ice velocity point (cp_ice_msh)                 !       |
1018         ! (C-grid ==>(U,V) or B-grid ==> I or F)                 j      |   T   U
1019         !                                                               |       |
1020         !                                                   j    j-1   -I-------|
1021         !                                               (for I)         |       |
1022         !                                                              i-1  i   i
1023         !                                                               i      i+1 (for I)
1024         SELECT CASE ( cp_ice_msh )
1025            !
1026         CASE( 'I' )                                         ! B-grid ==> I
1027            SELECT CASE ( srcv(jpr_itx1)%clgrid )
1028            CASE( 'U' )
1029               DO jj = 2, jpjm1                                   ! (U,V) ==> I
1030                  DO ji = 2, jpim1   ! NO vector opt.
1031                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji-1,jj  ,1) + frcv(jpr_itx1)%z3(ji-1,jj-1,1) )
1032                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji  ,jj-1,1) + frcv(jpr_ity1)%z3(ji-1,jj-1,1) )
1033                  END DO
1034               END DO
1035            CASE( 'F' )
1036               DO jj = 2, jpjm1                                   ! F ==> I
1037                  DO ji = 2, jpim1   ! NO vector opt.
1038                     p_taui(ji,jj) = frcv(jpr_itx1)%z3(ji-1,jj-1,1)
1039                     p_tauj(ji,jj) = frcv(jpr_ity1)%z3(ji-1,jj-1,1)
1040                  END DO
1041               END DO
1042            CASE( 'T' )
1043               DO jj = 2, jpjm1                                   ! T ==> I
1044                  DO ji = 2, jpim1   ! NO vector opt.
1045                     p_taui(ji,jj) = 0.25 * ( frcv(jpr_itx1)%z3(ji,jj  ,1) + frcv(jpr_itx1)%z3(ji-1,jj  ,1)   &
1046                        &                   + frcv(jpr_itx1)%z3(ji,jj-1,1) + frcv(jpr_itx1)%z3(ji-1,jj-1,1) ) 
1047                     p_tauj(ji,jj) = 0.25 * ( frcv(jpr_ity1)%z3(ji,jj  ,1) + frcv(jpr_ity1)%z3(ji-1,jj  ,1)   &
1048                        &                   + frcv(jpr_oty1)%z3(ji,jj-1,1) + frcv(jpr_ity1)%z3(ji-1,jj-1,1) )
1049                  END DO
1050               END DO
1051            CASE( 'I' )
1052               p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! I ==> I
1053               p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1054            END SELECT
1055            IF( srcv(jpr_itx1)%clgrid /= 'I' ) THEN
1056               CALL lbc_lnk( p_taui, 'I',  -1. )   ;   CALL lbc_lnk( p_tauj, 'I',  -1. )
1057            ENDIF
1058            !
1059         CASE( 'F' )                                         ! B-grid ==> F
1060            SELECT CASE ( srcv(jpr_itx1)%clgrid )
1061            CASE( 'U' )
1062               DO jj = 2, jpjm1                                   ! (U,V) ==> F
1063                  DO ji = fs_2, fs_jpim1   ! vector opt.
1064                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji,jj,1) + frcv(jpr_itx1)%z3(ji  ,jj+1,1) )
1065                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji,jj,1) + frcv(jpr_ity1)%z3(ji+1,jj  ,1) )
1066                  END DO
1067               END DO
1068            CASE( 'I' )
1069               DO jj = 2, jpjm1                                   ! I ==> F
1070                  DO ji = 2, jpim1   ! NO vector opt.
1071                     p_taui(ji,jj) = frcv(jpr_itx1)%z3(ji+1,jj+1,1)
1072                     p_tauj(ji,jj) = frcv(jpr_ity1)%z3(ji+1,jj+1,1)
1073                  END DO
1074               END DO
1075            CASE( 'T' )
1076               DO jj = 2, jpjm1                                   ! T ==> F
1077                  DO ji = 2, jpim1   ! NO vector opt.
1078                     p_taui(ji,jj) = 0.25 * ( frcv(jpr_itx1)%z3(ji,jj  ,1) + frcv(jpr_itx1)%z3(ji+1,jj  ,1)   &
1079                        &                   + frcv(jpr_itx1)%z3(ji,jj+1,1) + frcv(jpr_itx1)%z3(ji+1,jj+1,1) ) 
1080                     p_tauj(ji,jj) = 0.25 * ( frcv(jpr_ity1)%z3(ji,jj  ,1) + frcv(jpr_ity1)%z3(ji+1,jj  ,1)   &
1081                        &                   + frcv(jpr_ity1)%z3(ji,jj+1,1) + frcv(jpr_ity1)%z3(ji+1,jj+1,1) )
1082                  END DO
1083               END DO
1084            CASE( 'F' )
1085               p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! F ==> F
1086               p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1087            END SELECT
1088            IF( srcv(jpr_itx1)%clgrid /= 'F' ) THEN
1089               CALL lbc_lnk( p_taui, 'F',  -1. )   ;   CALL lbc_lnk( p_tauj, 'F',  -1. )
1090            ENDIF
1091            !
1092         CASE( 'C' )                                         ! C-grid ==> U,V
1093            SELECT CASE ( srcv(jpr_itx1)%clgrid )
1094            CASE( 'U' )
1095               p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! (U,V) ==> (U,V)
1096               p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1097            CASE( 'F' )
1098               DO jj = 2, jpjm1                                   ! F ==> (U,V)
1099                  DO ji = fs_2, fs_jpim1   ! vector opt.
1100                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji,jj,1) + frcv(jpr_itx1)%z3(ji  ,jj-1,1) )
1101                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(jj,jj,1) + frcv(jpr_ity1)%z3(ji-1,jj  ,1) )
1102                  END DO
1103               END DO
1104            CASE( 'T' )
1105               DO jj = 2, jpjm1                                   ! T ==> (U,V)
1106                  DO ji = fs_2, fs_jpim1   ! vector opt.
1107                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj  ,1) + frcv(jpr_itx1)%z3(ji,jj,1) )
1108                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji  ,jj+1,1) + frcv(jpr_ity1)%z3(ji,jj,1) )
1109                  END DO
1110               END DO
1111            CASE( 'I' )
1112               DO jj = 2, jpjm1                                   ! I ==> (U,V)
1113                  DO ji = 2, jpim1   ! NO vector opt.
1114                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj+1,1) + frcv(jpr_itx1)%z3(ji+1,jj  ,1) )
1115                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji+1,jj+1,1) + frcv(jpr_ity1)%z3(ji  ,jj+1,1) )
1116                  END DO
1117               END DO
1118            END SELECT
1119            IF( srcv(jpr_itx1)%clgrid /= 'U' ) THEN
1120               CALL lbc_lnk( p_taui, 'U',  -1. )   ;   CALL lbc_lnk( p_tauj, 'V',  -1. )
1121            ENDIF
1122         END SELECT
1123
1124      ENDIF
1125      !   
1126      CALL wrk_dealloc( jpi,jpj, ztx, zty )
1127      !
1128      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_ice_tau')
1129      !
1130   END SUBROUTINE sbc_cpl_ice_tau
1131   
1132
1133   SUBROUTINE sbc_cpl_ice_flx( p_frld  , palbi   , psst    , pist    )
1134      !!----------------------------------------------------------------------
1135      !!             ***  ROUTINE sbc_cpl_ice_flx  ***
1136      !!
1137      !! ** Purpose :   provide the heat and freshwater fluxes of the
1138      !!              ocean-ice system.
1139      !!
1140      !! ** Method  :   transform the fields received from the atmosphere into
1141      !!             surface heat and fresh water boundary condition for the
1142      !!             ice-ocean system. The following fields are provided:
1143      !!              * total non solar, solar and freshwater fluxes (qns_tot,
1144      !!             qsr_tot and emp_tot) (total means weighted ice-ocean flux)
1145      !!             NB: emp_tot include runoffs and calving.
1146      !!              * fluxes over ice (qns_ice, qsr_ice, emp_ice) where
1147      !!             emp_ice = sublimation - solid precipitation as liquid
1148      !!             precipitation are re-routed directly to the ocean and
1149      !!             runoffs and calving directly enter the ocean.
1150      !!              * solid precipitation (sprecip), used to add to qns_tot
1151      !!             the heat lost associated to melting solid precipitation
1152      !!             over the ocean fraction.
1153      !!       ===>> CAUTION here this changes the net heat flux received from
1154      !!             the atmosphere
1155      !!
1156      !!                  - the fluxes have been separated from the stress as
1157      !!                 (a) they are updated at each ice time step compare to
1158      !!                 an update at each coupled time step for the stress, and
1159      !!                 (b) the conservative computation of the fluxes over the
1160      !!                 sea-ice area requires the knowledge of the ice fraction
1161      !!                 after the ice advection and before the ice thermodynamics,
1162      !!                 so that the stress is updated before the ice dynamics
1163      !!                 while the fluxes are updated after it.
1164      !!
1165      !! ** Action  :   update at each nf_ice time step:
1166      !!                   qns_tot, qsr_tot  non-solar and solar total heat fluxes
1167      !!                   qns_ice, qsr_ice  non-solar and solar heat fluxes over the ice
1168      !!                   emp_tot            total evaporation - precipitation(liquid and solid) (-runoff)(-calving)
1169      !!                   emp_ice            ice sublimation - solid precipitation over the ice
1170      !!                   dqns_ice           d(non-solar heat flux)/d(Temperature) over the ice
1171      !!                   sprecip             solid precipitation over the ocean 
1172      !!----------------------------------------------------------------------
1173      REAL(wp), INTENT(in   ), DIMENSION(:,:)   ::   p_frld     ! lead fraction                [0 to 1]
1174      ! optional arguments, used only in 'mixed oce-ice' case
1175      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   palbi   ! all skies ice albedo
1176      REAL(wp), INTENT(in   ), DIMENSION(:,:  ), OPTIONAL ::   psst    ! sea surface temperature     [Celsius]
1177      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   pist    ! ice surface temperature     [Kelvin]
1178      !
1179      INTEGER ::   jl   ! dummy loop index
1180      REAL(wp), POINTER, DIMENSION(:,:) ::   zcptn, ztmp, zicefr
1181      !!----------------------------------------------------------------------
1182      !
1183      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_flx')
1184      !
1185      CALL wrk_alloc( jpi,jpj, zcptn, ztmp, zicefr )
1186
1187      zicefr(:,:) = 1.- p_frld(:,:)
1188      zcptn(:,:) = rcp * sst_m(:,:)
1189      !
1190      !                                                      ! ========================= !
1191      !                                                      !    freshwater budget      !   (emp)
1192      !                                                      ! ========================= !
1193      !
1194      !                                                           ! total Precipitations - total Evaporation (emp_tot)
1195      !                                                           ! solid precipitation  - sublimation       (emp_ice)
1196      !                                                           ! solid Precipitation                      (sprecip)
1197      SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
1198      CASE( 'conservative'  )   ! received fields: jpr_rain, jpr_snow, jpr_ievp, jpr_tevp
1199         sprecip(:,:) = frcv(jpr_snow)%z3(:,:,1)                 ! May need to ensure positive here
1200         tprecip(:,:) = frcv(jpr_rain)%z3(:,:,1) + sprecip (:,:) ! May need to ensure positive here
1201         emp_tot(:,:) = frcv(jpr_tevp)%z3(:,:,1) - tprecip(:,:)
1202#if defined key_cice
1203         IF ( TRIM(sn_rcv_emp%clcat) == 'yes' ) THEN
1204            ! emp_ice is the sum of frcv(jpr_ievp)%z3(:,:,1) over all layers - snow
1205            emp_ice(:,:) = - frcv(jpr_snow)%z3(:,:,1)
1206            DO jl=1,jpl
1207               emp_ice(:,:   ) = emp_ice(:,:) + frcv(jpr_ievp)%z3(:,:,jl)
1208            ENDDO
1209            ! latent heat coupled for each category in CICE
1210            qla_ice(:,:,1:jpl) = - frcv(jpr_ievp)%z3(:,:,1:jpl) * lsub
1211         ELSE
1212            ! If CICE has multicategories it still expects coupling fields for
1213            ! each even if we treat as a single field
1214            ! The latent heat flux is split between the ice categories according
1215            ! to the fraction of the ice in each category
1216            emp_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1)
1217            WHERE ( zicefr(:,:) /= 0._wp ) 
1218               ztmp(:,:) = 1./zicefr(:,:)
1219            ELSEWHERE 
1220               ztmp(:,:) = 0.e0
1221            END WHERE 
1222            DO jl=1,jpl
1223               qla_ice(:,:,jl) = - a_i(:,:,jl) * ztmp(:,:) * frcv(jpr_ievp)%z3(:,:,1) * lsub 
1224            END DO
1225            WHERE ( zicefr(:,:) == 0._wp )  qla_ice(:,:,1) = -frcv(jpr_ievp)%z3(:,:,1) * lsub 
1226         ENDIF
1227#else
1228         emp_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1)
1229#endif
1230            CALL iom_put( 'rain'         , frcv(jpr_rain)%z3(:,:,1)              )   ! liquid precipitation
1231         IF( iom_use('hflx_rain_cea') )   &
1232            CALL iom_put( 'hflx_rain_cea', frcv(jpr_rain)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from liq. precip.
1233         IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') )   &
1234            ztmp(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:)
1235         IF( iom_use('evap_ao_cea'  ) )   &
1236            CALL iom_put( 'evap_ao_cea'  , ztmp                   )   ! ice-free oce evap (cell average)
1237         IF( iom_use('hflx_evap_cea') )   &
1238            CALL iom_put( 'hflx_evap_cea', ztmp(:,:) * zcptn(:,:) )   ! heat flux from from evap (cell average)
1239      CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
1240         emp_tot(:,:) = p_frld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + zicefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
1241         emp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1)
1242         sprecip(:,:) = - frcv(jpr_semp)%z3(:,:,1) + frcv(jpr_ievp)%z3(:,:,1)
1243      END SELECT
1244
1245         CALL iom_put( 'snowpre'    , sprecip                                )   ! Snow
1246      IF( iom_use('snow_ao_cea') )   &
1247         CALL iom_put( 'snow_ao_cea', sprecip(:,:) * p_frld(:,:)             )   ! Snow        over ice-free ocean  (cell average)
1248      IF( iom_use('snow_ai_cea') )   &
1249         CALL iom_put( 'snow_ai_cea', sprecip(:,:) * zicefr(:,:)             )   ! Snow        over sea-ice         (cell average)
1250      IF( iom_use('subl_ai_cea') )   &
1251         CALL iom_put( 'subl_ai_cea', frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) )   ! Sublimation over sea-ice         (cell average)
1252      !   
1253      !                                                           ! runoffs and calving (put in emp_tot)
1254      IF( srcv(jpr_rnf)%laction ) THEN
1255         emp_tot(:,:) = emp_tot(:,:) - frcv(jpr_rnf)%z3(:,:,1)
1256            CALL iom_put( 'runoffs'      , frcv(jpr_rnf)%z3(:,:,1)              )   ! rivers
1257         IF( iom_use('hflx_rnf_cea') )   &
1258            CALL iom_put( 'hflx_rnf_cea' , frcv(jpr_rnf)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from rivers
1259      ENDIF
1260      IF( srcv(jpr_cal)%laction ) THEN
1261         emp_tot(:,:) = emp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
1262         CALL iom_put( 'calving', frcv(jpr_cal)%z3(:,:,1) )
1263      ENDIF
1264      !
1265!!gm :  this seems to be internal cooking, not sure to need that in a generic interface
1266!!gm                                       at least should be optional...
1267!!       ! remove negative runoff                            ! sum over the global domain
1268!!       zcumulpos = SUM( MAX( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * e1t(:,:) * e2t(:,:) * tmask_i(:,:) )
1269!!       zcumulneg = SUM( MIN( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * e1t(:,:) * e2t(:,:) * tmask_i(:,:) )
1270!!       IF( lk_mpp )   CALL mpp_sum( zcumulpos )
1271!!       IF( lk_mpp )   CALL mpp_sum( zcumulneg )
1272!!       IF( zcumulpos /= 0. ) THEN                          ! distribute negative runoff on positive runoff grid points
1273!!          zcumulneg = 1.e0 + zcumulneg / zcumulpos
1274!!          frcv(jpr_rnf)%z3(:,:,1) = MAX( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * zcumulneg
1275!!       ENDIF     
1276!!       emp_tot(:,:) = emp_tot(:,:) - frcv(jpr_rnf)%z3(:,:,1)   ! add runoff to e-p
1277!!
1278!!gm  end of internal cooking
1279
1280      !                                                      ! ========================= !
1281      SELECT CASE( TRIM( sn_rcv_qns%cldes ) )                !   non solar heat fluxes   !   (qns)
1282      !                                                      ! ========================= !
1283      CASE( 'oce only' )                                     ! the required field is directly provided
1284         qns_tot(:,:  ) = frcv(jpr_qnsoce)%z3(:,:,1)
1285      CASE( 'conservative' )                                      ! the required fields are directly provided
1286         qns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
1287         IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
1288            qns_ice(:,:,1:jpl) = frcv(jpr_qnsice)%z3(:,:,1:jpl)
1289         ELSE
1290            ! Set all category values equal for the moment
1291            DO jl=1,jpl
1292               qns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
1293            ENDDO
1294         ENDIF
1295      CASE( 'oce and ice' )       ! the total flux is computed from ocean and ice fluxes
1296         qns_tot(:,:  ) =  p_frld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1)
1297         IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
1298            DO jl=1,jpl
1299               qns_tot(:,:   ) = qns_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qnsice)%z3(:,:,jl)   
1300               qns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,jl)
1301            ENDDO
1302         ELSE
1303            DO jl=1,jpl
1304               qns_tot(:,:   ) = qns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
1305               qns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
1306            ENDDO
1307         ENDIF
1308      CASE( 'mixed oce-ice' )     ! the ice flux is cumputed from the total flux, the SST and ice informations
1309! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
1310         qns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
1311         qns_ice(:,:,1) = frcv(jpr_qnsmix)%z3(:,:,1)    &
1312            &            + frcv(jpr_dqnsdt)%z3(:,:,1) * ( pist(:,:,1) - ( (rt0 + psst(:,:  ) ) * p_frld(:,:)   &
1313            &                                                   +          pist(:,:,1)   * zicefr(:,:) ) )
1314      END SELECT
1315      ztmp(:,:) = p_frld(:,:) * sprecip(:,:) * lfus
1316      qns_tot(:,:) = qns_tot(:,:)                         &            ! qns_tot update over free ocean with:
1317         &          - ztmp(:,:)                           &            ! remove the latent heat flux of solid precip. melting
1318         &          - (  emp_tot(:,:)                     &            ! remove the heat content of mass flux (assumed to be at SST)
1319         &             - emp_ice(:,:) * zicefr(:,:)  ) * zcptn(:,:) 
1320      IF( iom_use('hflx_snow_cea') )   &
1321         CALL iom_put( 'hflx_snow_cea', ztmp + sprecip(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average)
1322!!gm
1323!!    currently it is taken into account in leads budget but not in the qns_tot, and thus not in
1324!!    the flux that enter the ocean....
1325!!    moreover 1 - it is not diagnose anywhere....
1326!!             2 - it is unclear for me whether this heat lost is taken into account in the atmosphere or not...
1327!!
1328!! similar job should be done for snow and precipitation temperature
1329      !                                     
1330      IF( srcv(jpr_cal)%laction ) THEN                            ! Iceberg melting
1331         ztmp(:,:) = frcv(jpr_cal)%z3(:,:,1) * lfus               ! add the latent heat of iceberg melting
1332         qns_tot(:,:) = qns_tot(:,:) - ztmp(:,:)
1333         IF( iom_use('hflx_cal_cea') )   &
1334            CALL iom_put( 'hflx_cal_cea', ztmp + frcv(jpr_cal)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from calving
1335      ENDIF
1336
1337      !                                                      ! ========================= !
1338      SELECT CASE( TRIM( sn_rcv_qsr%cldes ) )                !      solar heat fluxes    !   (qsr)
1339      !                                                      ! ========================= !
1340      CASE( 'oce only' )
1341         qsr_tot(:,:  ) = MAX( 0._wp , frcv(jpr_qsroce)%z3(:,:,1) )
1342      CASE( 'conservative' )
1343         qsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
1344         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
1345            qsr_ice(:,:,1:jpl) = frcv(jpr_qsrice)%z3(:,:,1:jpl)
1346         ELSE
1347            ! Set all category values equal for the moment
1348            DO jl=1,jpl
1349               qsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
1350            ENDDO
1351         ENDIF
1352         qsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
1353         qsr_ice(:,:,1) = frcv(jpr_qsrice)%z3(:,:,1)
1354      CASE( 'oce and ice' )
1355         qsr_tot(:,:  ) =  p_frld(:,:) * frcv(jpr_qsroce)%z3(:,:,1)
1356         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
1357            DO jl=1,jpl
1358               qsr_tot(:,:   ) = qsr_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qsrice)%z3(:,:,jl)   
1359               qsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,jl)
1360            ENDDO
1361         ELSE
1362            DO jl=1,jpl
1363               qsr_tot(:,:   ) = qsr_tot(:,:) + zicefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
1364               qsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
1365            ENDDO
1366         ENDIF
1367      CASE( 'mixed oce-ice' )
1368         qsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
1369! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
1370!       Create solar heat flux over ice using incoming solar heat flux and albedos
1371!       ( see OASIS3 user guide, 5th edition, p39 )
1372         qsr_ice(:,:,1) = frcv(jpr_qsrmix)%z3(:,:,1) * ( 1.- palbi(:,:,1) )   &
1373            &            / (  1.- ( albedo_oce_mix(:,:  ) * p_frld(:,:)       &
1374            &                     + palbi         (:,:,1) * zicefr(:,:) ) )
1375      END SELECT
1376      IF( ln_dm2dc ) THEN   ! modify qsr to include the diurnal cycle
1377         qsr_tot(:,:  ) = sbc_dcy( qsr_tot(:,:  ) )
1378         DO jl=1,jpl
1379            qsr_ice(:,:,jl) = sbc_dcy( qsr_ice(:,:,jl) )
1380         ENDDO
1381      ENDIF
1382
1383      !                                                      ! ========================= !
1384      SELECT CASE( TRIM( sn_rcv_dqnsdt%cldes ) )             !          d(qns)/dt        !
1385      !                                                      ! ========================= !
1386      CASE ('coupled')
1387         IF ( TRIM(sn_rcv_dqnsdt%clcat) == 'yes' ) THEN
1388            dqns_ice(:,:,1:jpl) = frcv(jpr_dqnsdt)%z3(:,:,1:jpl)
1389         ELSE
1390            ! Set all category values equal for the moment
1391            DO jl=1,jpl
1392               dqns_ice(:,:,jl) = frcv(jpr_dqnsdt)%z3(:,:,1)
1393            ENDDO
1394         ENDIF
1395      END SELECT
1396
1397      !                                                      ! ========================= !
1398      SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) )             !    topmelt and botmelt    !
1399      !                                                      ! ========================= !
1400      CASE ('coupled')
1401         topmelt(:,:,:)=frcv(jpr_topm)%z3(:,:,:)
1402         botmelt(:,:,:)=frcv(jpr_botm)%z3(:,:,:)
1403      END SELECT
1404
1405      ! Surface transimission parameter io (Maykut Untersteiner , 1971 ; Ebert and Curry, 1993 )
1406      ! Used for LIM2 and LIM3
1407      ! Coupled case: since cloud cover is not received from atmosphere
1408      !               ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
1409      fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
1410      fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
1411
1412      CALL wrk_dealloc( jpi,jpj, zcptn, ztmp, zicefr )
1413      !
1414      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_ice_flx')
1415      !
1416   END SUBROUTINE sbc_cpl_ice_flx
1417   
1418   
1419   SUBROUTINE sbc_cpl_snd( kt )
1420      !!----------------------------------------------------------------------
1421      !!             ***  ROUTINE sbc_cpl_snd  ***
1422      !!
1423      !! ** Purpose :   provide the ocean-ice informations to the atmosphere
1424      !!
1425      !! ** Method  :   send to the atmosphere through a call to cpl_snd
1426      !!              all the needed fields (as defined in sbc_cpl_init)
1427      !!----------------------------------------------------------------------
1428      INTEGER, INTENT(in) ::   kt
1429      !
1430      INTEGER ::   ji, jj, jl   ! dummy loop indices
1431      INTEGER ::   isec, info   ! local integer
1432      REAL(wp), POINTER, DIMENSION(:,:)   ::   zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1
1433      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmp3, ztmp4   
1434      !!----------------------------------------------------------------------
1435      !
1436      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_snd')
1437      !
1438      CALL wrk_alloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
1439      CALL wrk_alloc( jpi,jpj,jpl, ztmp3, ztmp4 )
1440
1441      isec = ( kt - nit000 ) * NINT(rdttra(1))        ! date of exchanges
1442
1443      zfr_l(:,:) = 1.- fr_i(:,:)
1444      !                                                      ! ------------------------- !
1445      !                                                      !    Surface temperature    !   in Kelvin
1446      !                                                      ! ------------------------- !
1447      IF( ssnd(jps_toce)%laction .OR. ssnd(jps_tice)%laction .OR. ssnd(jps_tmix)%laction ) THEN
1448         SELECT CASE( sn_snd_temp%cldes)
1449         CASE( 'oce only'             )   ;   ztmp1(:,:) =   tsn(:,:,1,jp_tem) + rt0
1450         CASE( 'weighted oce and ice' )   ;   ztmp1(:,:) = ( tsn(:,:,1,jp_tem) + rt0 ) * zfr_l(:,:)   
1451            SELECT CASE( sn_snd_temp%clcat )
1452            CASE( 'yes' )   
1453               ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
1454            CASE( 'no' )
1455               ztmp3(:,:,:) = 0.0
1456               DO jl=1,jpl
1457                  ztmp3(:,:,1) = ztmp3(:,:,1) + tn_ice(:,:,jl) * a_i(:,:,jl)
1458               ENDDO
1459            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
1460            END SELECT
1461         CASE( 'mixed oce-ice'        )   
1462            ztmp1(:,:) = ( tsn(:,:,1,jp_tem) + rt0 ) * zfr_l(:,:) 
1463            DO jl=1,jpl
1464               ztmp1(:,:) = ztmp1(:,:) + tn_ice(:,:,jl) * a_i(:,:,jl)
1465            ENDDO
1466         CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%cldes' )
1467         END SELECT
1468         IF( ssnd(jps_toce)%laction )   CALL cpl_snd( jps_toce, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
1469         IF( ssnd(jps_tice)%laction )   CALL cpl_snd( jps_tice, isec, ztmp3, info )
1470         IF( ssnd(jps_tmix)%laction )   CALL cpl_snd( jps_tmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
1471      ENDIF
1472      !                                                      ! ------------------------- !
1473      !                                                      !           Albedo          !
1474      !                                                      ! ------------------------- !
1475      IF( ssnd(jps_albice)%laction ) THEN                         ! ice
1476         ztmp3(:,:,1:jpl) = alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
1477         CALL cpl_snd( jps_albice, isec, ztmp3, info )
1478      ENDIF
1479      IF( ssnd(jps_albmix)%laction ) THEN                         ! mixed ice-ocean
1480         ztmp1(:,:) = albedo_oce_mix(:,:) * zfr_l(:,:)
1481         DO jl=1,jpl
1482            ztmp1(:,:) = ztmp1(:,:) + alb_ice(:,:,jl) * a_i(:,:,jl)
1483         ENDDO
1484         CALL cpl_snd( jps_albmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
1485      ENDIF
1486      !                                                      ! ------------------------- !
1487      !                                                      !  Ice fraction & Thickness !
1488      !                                                      ! ------------------------- !
1489      ! Send ice fraction field
1490      IF( ssnd(jps_fice)%laction ) THEN
1491         SELECT CASE( sn_snd_thick%clcat )
1492         CASE( 'yes' )   ;   ztmp3(:,:,1:jpl) =  a_i(:,:,1:jpl)
1493         CASE( 'no'  )   ;   ztmp3(:,:,1    ) = fr_i(:,:      )
1494         CASE default    ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
1495         END SELECT
1496         CALL cpl_snd( jps_fice, isec, ztmp3, info )
1497      ENDIF
1498
1499      ! Send ice and snow thickness field
1500      IF( ssnd(jps_hice)%laction .OR. ssnd(jps_hsnw)%laction ) THEN
1501         SELECT CASE( sn_snd_thick%cldes)
1502         CASE( 'none'                  )       ! nothing to do
1503         CASE( 'weighted ice and snow' )   
1504            SELECT CASE( sn_snd_thick%clcat )
1505            CASE( 'yes' )   
1506               ztmp3(:,:,1:jpl) =  ht_i(:,:,1:jpl) * a_i(:,:,1:jpl)
1507               ztmp4(:,:,1:jpl) =  ht_s(:,:,1:jpl) * a_i(:,:,1:jpl)
1508            CASE( 'no' )
1509               ztmp3(:,:,:) = 0.0   ;  ztmp4(:,:,:) = 0.0
1510               DO jl=1,jpl
1511                  ztmp3(:,:,1) = ztmp3(:,:,1) + ht_i(:,:,jl) * a_i(:,:,jl)
1512                  ztmp4(:,:,1) = ztmp4(:,:,1) + ht_s(:,:,jl) * a_i(:,:,jl)
1513               ENDDO
1514            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
1515            END SELECT
1516         CASE( 'ice and snow'         )   
1517            ztmp3(:,:,1:jpl) = ht_i(:,:,1:jpl)
1518            ztmp4(:,:,1:jpl) = ht_s(:,:,1:jpl)
1519         CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%cldes' )
1520         END SELECT
1521         IF( ssnd(jps_hice)%laction )   CALL cpl_snd( jps_hice, isec, ztmp3, info )
1522         IF( ssnd(jps_hsnw)%laction )   CALL cpl_snd( jps_hsnw, isec, ztmp4, info )
1523      ENDIF
1524      !
1525      ! Send meltpond fields
1526      IF( ssnd(jps_a_p)%laction .OR. ssnd(jps_ht_p)%laction ) THEN
1527         SELECT CASE( sn_snd_mpnd%cldes) 
1528         CASE( 'weighted ice' ) 
1529            SELECT CASE( sn_snd_mpnd%clcat ) 
1530            CASE( 'yes' ) 
1531               ztmp3(:,:,1:jpl) =  a_p(:,:,1:jpl) * a_i(:,:,1:jpl) 
1532               ztmp4(:,:,1:jpl) =  ht_p(:,:,1:jpl) * a_i(:,:,1:jpl) 
1533            CASE( 'no' ) 
1534               ztmp3(:,:,:) = 0.0 
1535               ztmp4(:,:,:) = 0.0 
1536               DO jl=1,jpl 
1537                 ztmp3(:,:,1) = ztmp3(:,:,1) + a_p(:,:,jpl) * a_i(:,:,jpl) 
1538                 ztmp4(:,:,1) = ztmp4(:,:,1) + ht_p(:,:,jpl) * a_i(:,:,jpl) 
1539               ENDDO 
1540            CASE default    ;   CALL ctl_stop( 'sbc_cpl_mpd: wrong definition of sn_snd_mpnd%clcat' ) 
1541            END SELECT
1542         CASE( 'ice only' )   
1543            ztmp3(:,:,1:jpl) = a_p(:,:,1:jpl) 
1544            ztmp4(:,:,1:jpl) = ht_p(:,:,1:jpl) 
1545         END SELECT
1546         IF( ssnd(jps_a_p)%laction )   CALL cpl_snd( jps_a_p, isec, ztmp3, info )   
1547         IF( ssnd(jps_ht_p)%laction )   CALL cpl_snd( jps_ht_p, isec, ztmp4, info )   
1548      ENDIF
1549      !
1550      !
1551#if defined key_cpl_carbon_cycle
1552      !                                                      ! ------------------------- !
1553      !                                                      !  CO2 flux from PISCES     !
1554      !                                                      ! ------------------------- !
1555      IF( ssnd(jps_co2)%laction )   CALL cpl_snd( jps_co2, isec, RESHAPE ( oce_co2, (/jpi,jpj,1/) ) , info )
1556      !
1557#endif
1558      !                                                      ! ------------------------- !
1559      IF( ssnd(jps_ocx1)%laction ) THEN                      !      Surface current      !
1560         !                                                   ! ------------------------- !
1561         !   
1562         !                                                  j+1   j     -----V---F
1563         ! surface velocity always sent from T point                     !       |
1564         !                                                        j      |   T   U
1565         !                                                               |       |
1566         !                                                   j    j-1   -I-------|
1567         !                                               (for I)         |       |
1568         !                                                              i-1  i   i
1569         !                                                               i      i+1 (for I)
1570         SELECT CASE( TRIM( sn_snd_crt%cldes ) )
1571         CASE( 'oce only'             )      ! C-grid ==> T
1572            DO jj = 2, jpjm1
1573               DO ji = fs_2, fs_jpim1   ! vector opt.
1574                  zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj  ,1) )
1575                  zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji  ,jj-1,1) ) 
1576               END DO
1577            END DO
1578         CASE( 'weighted oce and ice' )   
1579            SELECT CASE ( cp_ice_msh )
1580            CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
1581               DO jj = 2, jpjm1
1582                  DO ji = fs_2, fs_jpim1   ! vector opt.
1583                     zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj) 
1584                     zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)
1585                     zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
1586                     zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
1587                  END DO
1588               END DO
1589            CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
1590               DO jj = 2, jpjm1
1591                  DO ji = 2, jpim1   ! NO vector opt.
1592                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj) 
1593                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj) 
1594                     zitx1(ji,jj) = 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &
1595                        &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1596                     zity1(ji,jj) = 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &
1597                        &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1598                  END DO
1599               END DO
1600            CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
1601               DO jj = 2, jpjm1
1602                  DO ji = 2, jpim1   ! NO vector opt.
1603                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj) 
1604                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj) 
1605                     zitx1(ji,jj) = 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &
1606                        &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1607                     zity1(ji,jj) = 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &
1608                        &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1609                  END DO
1610               END DO
1611            END SELECT
1612            CALL lbc_lnk( zitx1, 'T', -1. )   ;   CALL lbc_lnk( zity1, 'T', -1. )
1613         CASE( 'mixed oce-ice'        )
1614            SELECT CASE ( cp_ice_msh )
1615            CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
1616               DO jj = 2, jpjm1
1617                  DO ji = fs_2, fs_jpim1   ! vector opt.
1618                     zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)   &
1619                        &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
1620                     zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)   &
1621                        &         + 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
1622                  END DO
1623               END DO
1624            CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
1625               DO jj = 2, jpjm1
1626                  DO ji = 2, jpim1   ! NO vector opt.
1627                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
1628                        &         + 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &
1629                        &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1630                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
1631                        &         + 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &
1632                        &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1633                  END DO
1634               END DO
1635            CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
1636               DO jj = 2, jpjm1
1637                  DO ji = 2, jpim1   ! NO vector opt.
1638                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
1639                        &         + 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &
1640                        &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1641                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
1642                        &         + 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &
1643                        &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1644                  END DO
1645               END DO
1646            END SELECT
1647         END SELECT
1648         CALL lbc_lnk( zotx1, ssnd(jps_ocx1)%clgrid, -1. )   ;   CALL lbc_lnk( zoty1, ssnd(jps_ocy1)%clgrid, -1. )
1649         !
1650         !
1651         IF( TRIM( sn_snd_crt%clvor ) == 'eastward-northward' ) THEN             ! Rotation of the components
1652            !                                                                     ! Ocean component
1653            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->e', ztmp1 )       ! 1st component
1654            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->n', ztmp2 )       ! 2nd component
1655            zotx1(:,:) = ztmp1(:,:)                                                   ! overwrite the components
1656            zoty1(:,:) = ztmp2(:,:)
1657            IF( ssnd(jps_ivx1)%laction ) THEN                                     ! Ice component
1658               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 )    ! 1st component
1659               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 )    ! 2nd component
1660               zitx1(:,:) = ztmp1(:,:)                                                ! overwrite the components
1661               zity1(:,:) = ztmp2(:,:)
1662            ENDIF
1663         ENDIF
1664         !
1665         ! spherical coordinates to cartesian -> 2 components to 3 components
1666         IF( TRIM( sn_snd_crt%clvref ) == 'cartesian' ) THEN
1667            ztmp1(:,:) = zotx1(:,:)                     ! ocean currents
1668            ztmp2(:,:) = zoty1(:,:)
1669            CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 )
1670            !
1671            IF( ssnd(jps_ivx1)%laction ) THEN           ! ice velocities
1672               ztmp1(:,:) = zitx1(:,:)
1673               ztmp1(:,:) = zity1(:,:)
1674               CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 )
1675            ENDIF
1676         ENDIF
1677         !
1678         IF( ssnd(jps_ocx1)%laction )   CALL cpl_snd( jps_ocx1, isec, RESHAPE ( zotx1, (/jpi,jpj,1/) ), info )   ! ocean x current 1st grid
1679         IF( ssnd(jps_ocy1)%laction )   CALL cpl_snd( jps_ocy1, isec, RESHAPE ( zoty1, (/jpi,jpj,1/) ), info )   ! ocean y current 1st grid
1680         IF( ssnd(jps_ocz1)%laction )   CALL cpl_snd( jps_ocz1, isec, RESHAPE ( zotz1, (/jpi,jpj,1/) ), info )   ! ocean z current 1st grid
1681         !
1682         IF( ssnd(jps_ivx1)%laction )   CALL cpl_snd( jps_ivx1, isec, RESHAPE ( zitx1, (/jpi,jpj,1/) ), info )   ! ice   x current 1st grid
1683         IF( ssnd(jps_ivy1)%laction )   CALL cpl_snd( jps_ivy1, isec, RESHAPE ( zity1, (/jpi,jpj,1/) ), info )   ! ice   y current 1st grid
1684         IF( ssnd(jps_ivz1)%laction )   CALL cpl_snd( jps_ivz1, isec, RESHAPE ( zitz1, (/jpi,jpj,1/) ), info )   ! ice   z current 1st grid
1685         !
1686      ENDIF
1687      !
1688      CALL wrk_dealloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
1689      CALL wrk_dealloc( jpi,jpj,jpl, ztmp3, ztmp4 )
1690      !
1691      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_snd')
1692      !
1693   END SUBROUTINE sbc_cpl_snd
1694   
1695   !!======================================================================
1696END MODULE sbccpl
Note: See TracBrowser for help on using the repository browser.