New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sbccpl.F90 in branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

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

Last change on this file since 5026 was 5026, checked in by timgraham, 9 years ago

First set of changes for multicategory coupling with CICE

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