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 @ 5119

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

Changes for coupling with meltponds

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