source: branches/UKMO/dev_r5107_hadgem3_mct/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 5284

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

First attempt to convert OASIS3-MCT branch from NEMO3.5 to NEMO3.6. Original branch was dev/frrh/vn3.5_beta_hadgem3_mct

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