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

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

source: branches/UKMO/dev_3841_sbc/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 4827

Last change on this file since 4827 was 4827, checked in by charris, 9 years ago

Some demonstration code changes.

  • Property svn:keywords set to Id
File size: 78.7 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#if defined key_oasis3 || defined key_oasis4
12   !!----------------------------------------------------------------------
13   !!   'key_oasis3' or 'key_oasis4'   Coupled Ocean/Atmosphere formulation
14   !!----------------------------------------------------------------------
15   !!   namsbc_cpl      : coupled formulation namlist
16   !!   sbc_cpl_init    : initialisation of the coupled exchanges
17   !!   sbc_cpl_rcv     : receive fields from the atmosphere over the ocean (ocean only)
18   !!                     receive stress from the atmosphere over the ocean (ocean-ice case)
19   !!   sbc_cpl_ice_tau : receive stress from the atmosphere over ice
20   !!   sbc_cpl_ice_flx : receive fluxes from the atmosphere over ice
21   !!   sbc_cpl_snd     : send     fields to the atmosphere
22   !!----------------------------------------------------------------------
23   USE dom_oce         ! ocean space and time domain
24   USE sbc_oce         ! Surface boundary condition: ocean fields
25   USE sbc_ice         ! Surface boundary condition: ice fields
26   USE sbcdcy          ! surface boundary condition: diurnal cycle
27   USE phycst          ! physical constants
28   USE fldread2, ONLY: fld_fill2           ! read input fields
29   USE fld_def
30   USE sbcget
31#if defined key_lim3
32   USE par_ice         ! ice parameters
33   USE ice             ! ice variables
34#endif
35#if defined key_lim2
36   USE par_ice_2       ! ice parameters
37   USE ice_2           ! ice variables
38#endif
39#if defined key_oasis3
40   USE cpl_oasis3      ! OASIS3 coupling
41#endif
42#if defined key_oasis4
43   USE cpl_oasis4      ! OASIS4 coupling
44#endif
45   USE geo2ocean       !
46   USE oce   , ONLY : tsn, un, vn
47   USE albedo          !
48   USE in_out_manager  ! I/O manager
49   USE iom             ! NetCDF library
50   USE lib_mpp         ! distribued memory computing library
51   USE wrk_nemo        ! work arrays
52   USE timing          ! Timing
53   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
54#if defined key_cpl_carbon_cycle
55   USE p4zflx, ONLY : oce_co2
56#endif
57   USE diaar5, ONLY :   lk_diaar5
58
59   IMPLICIT NONE
60   PRIVATE
61
62   PUBLIC   sbc_cpl_rcv        ! routine called by sbc_ice_lim(_2).F90
63   PUBLIC   sbc_cpl_snd        ! routine called by step.F90
64   PUBLIC   sbc_cpl_ice_tau    ! routine called by sbc_ice_lim(_2).F90
65   PUBLIC   sbc_cpl_ice_flx    ! routine called by sbc_ice_lim(_2).F90
66
67   INTEGER, PARAMETER ::   jps_fice   =  1            ! ice fraction
68   INTEGER, PARAMETER ::   jps_toce   =  2            ! ocean temperature
69   INTEGER, PARAMETER ::   jps_tice   =  3            ! ice   temperature
70   INTEGER, PARAMETER ::   jps_tmix   =  4            ! mixed temperature (ocean+ice)
71   INTEGER, PARAMETER ::   jps_albice =  5            ! ice   albedo
72   INTEGER, PARAMETER ::   jps_albmix =  6            ! mixed albedo
73   INTEGER, PARAMETER ::   jps_hice   =  7            ! ice  thickness
74   INTEGER, PARAMETER ::   jps_hsnw   =  8            ! snow thickness
75   INTEGER, PARAMETER ::   jps_ocx1   =  9            ! ocean current on grid 1
76   INTEGER, PARAMETER ::   jps_ocy1   = 10            !
77   INTEGER, PARAMETER ::   jps_ocz1   = 11            !
78   INTEGER, PARAMETER ::   jps_ivx1   = 12            ! ice   current on grid 1
79   INTEGER, PARAMETER ::   jps_ivy1   = 13            !
80   INTEGER, PARAMETER ::   jps_ivz1   = 14            !
81   INTEGER, PARAMETER ::   jps_co2    = 15
82   INTEGER, PARAMETER ::   jpsnd      = 15            ! total number of fields sended
83
84   !                                                         !!** namelist namsbc_cpl **
85   TYPE ::   FLD_C
86      CHARACTER(len = 32) ::   cldes                  ! desciption of the coupling strategy
87      CHARACTER(len = 32) ::   clcat                  ! multiple ice categories strategy
88      CHARACTER(len = 32) ::   clvref                 ! reference of vector ('spherical' or 'cartesian')
89      CHARACTER(len = 32) ::   clvor                  ! orientation of vector fields ('eastward-northward' or 'local grid')
90      CHARACTER(len = 32) ::   clvgrd                 ! grids on which is located the vector fields
91   END TYPE FLD_C
92
93
94   ! Send to the atmosphere                           !
95   TYPE(FLD_C) ::   sn_snd_temp, sn_snd_alb, sn_snd_thick, sn_snd_crt, sn_snd_co2                       
96   ! Received from the atmosphere                     !
97   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
98   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2                       
99   TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_otx1, sn_oty1, sn_otz1, sn_otx2, sn_oty2, sn_otz2,  sn_itx1, sn_ity1, sn_itz1, sn_itx2, sn_ity2, sn_itz2
100
101   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   albedo_oce_mix     ! ocean albedo sent to atmosphere (mix clear/overcast sky)
102
103#if ! defined key_lim2   &&   ! defined key_lim3
104   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice,fr1_i0,fr2_i0          ! jpi, jpj
105   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tn_ice, alb_ice, qns_ice, dqns_ice  ! (jpi,jpj,jpl)
106#endif
107
108#if ! defined key_cice  && ! defined key_lim2   &&   ! defined key_lim3
109   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   emp_ice
110   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qsr_ice
111#endif
112
113#if ! defined key_lim3   &&  ! defined key_cice
114   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  a_i
115#endif
116
117#if ! defined key_lim3
118   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  ht_i, ht_s
119#endif
120
121#if ! defined key_cice
122   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  topmelt, botmelt
123#endif
124
125   !! Substitution
126#  include "vectopt_loop_substitute.h90"
127   !!----------------------------------------------------------------------
128   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
129   !! $Id$
130   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
131   !!----------------------------------------------------------------------
132
133CONTAINS
134 
135   INTEGER FUNCTION sbc_cpl_alloc()
136      !!----------------------------------------------------------------------
137      !!             ***  FUNCTION sbc_cpl_alloc  ***
138      !!----------------------------------------------------------------------
139      INTEGER :: ierr(4),jn
140      !!----------------------------------------------------------------------
141      ierr(:) = 0
142      !
143      ALLOCATE( albedo_oce_mix(jpi,jpj),  STAT=ierr(1) )
144      !
145#if ! defined key_lim2 && ! defined key_lim3
146      ! quick patch to be able to run the coupled model without sea-ice...
147      ALLOCATE( u_ice(jpi,jpj) , fr1_i0(jpi,jpj) , tn_ice (jpi,jpj,1) ,     &
148                v_ice(jpi,jpj) , fr2_i0(jpi,jpj) , alb_ice(jpi,jpj,1),      &
149                emp_ice(jpi,jpj) , qns_ice(jpi,jpj,1) , dqns_ice(jpi,jpj,1) , STAT=ierr(2) )
150#endif
151
152#if ! defined key_lim3 && ! defined key_cice
153      ALLOCATE( a_i(jpi,jpj,jpl) , STAT=ierr(3) )
154#endif
155
156#if defined key_cice || defined key_lim2
157      ALLOCATE( ht_i(jpi,jpj,jpl) , ht_s(jpi,jpj,jpl) , STAT=ierr(4) )
158#endif
159      sbc_cpl_alloc = MAXVAL( ierr )
160      IF( lk_mpp            )   CALL mpp_sum ( sbc_cpl_alloc )
161      IF( sbc_cpl_alloc > 0 )   CALL ctl_warn('sbc_cpl_alloc: allocation of arrays failed')
162      !
163   END FUNCTION sbc_cpl_alloc
164
165
166   SUBROUTINE sbc_cpl_init( k_ice )     
167      !!----------------------------------------------------------------------
168      !!             ***  ROUTINE sbc_cpl_init  ***
169      !!
170      !! ** Purpose :   Initialisation of send and recieved information from
171      !!                the atmospheric component
172      !!
173      !! ** Method  : * Read namsbc_cpl namelist
174      !!              * define the receive interface
175      !!              * define the send    interface
176      !!              * initialise the OASIS coupler
177      !!----------------------------------------------------------------------
178      INTEGER, INTENT(in) ::   k_ice    ! ice management in the sbc (=0/1/2/3)
179      !!
180      INTEGER ::   jn       ! dummy loop index
181      REAL(wp), POINTER, DIMENSION(:,:) ::   zacs, zaos
182      !!
183      NAMELIST/namsbc_cpl/  sn_snd_temp, sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2,   &
184         &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr,   &
185         &                  sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf  , sn_rcv_cal   , sn_rcv_iceflx  , sn_rcv_co2
186      !!---------------------------------------------------------------------
187      !
188      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_init')
189      !
190      CALL wrk_alloc( jpi,jpj, zacs, zaos )
191
192      ! ================================ !
193      !      Namelist informations       !
194      ! ================================ !
195
196      ! default definitions
197      !                    !     description       !  multiple  !    vector   !      vector          ! vector !
198      !                    !                       ! categories !  reference  !    orientation       ! grids  !
199      ! send
200      sn_snd_temp   = FLD_C( 'weighted oce and ice',    'no'    ,     ''      ,         ''           ,   ''   ) 
201      sn_snd_alb    = FLD_C( 'weighted ice'        ,    'no'    ,     ''      ,         ''           ,   ''   ) 
202      sn_snd_thick  = FLD_C( 'none'                ,    'no'    ,     ''      ,         ''           ,   ''   ) 
203      sn_snd_crt    = FLD_C( 'none'                ,    'no'    , 'spherical' , 'eastward-northward' ,  'T'   )     
204      sn_snd_co2    = FLD_C( 'none'                ,    'no'    ,     ''      ,         ''           ,   ''   )     
205      ! receive
206      sn_rcv_w10m   = FLD_C( 'none'                ,    'no'    ,     ''      ,         ''          ,   ''    )
207      sn_rcv_taumod = FLD_C( 'coupled'             ,    'no'    ,     ''      ,         ''          ,   ''    )
208      sn_rcv_tau    = FLD_C( 'oce only'            ,    'no'    , 'cartesian' , 'eastward-northward',  'U,V'  ) 
209      sn_rcv_dqnsdt = FLD_C( 'coupled'             ,    'no'    ,     ''      ,         ''          ,   ''    )
210      sn_rcv_qsr    = FLD_C( 'oce and ice'         ,    'no'    ,     ''      ,         ''          ,   ''    )
211      sn_rcv_qns    = FLD_C( 'oce and ice'         ,    'no'    ,     ''      ,         ''          ,   ''    )
212      sn_rcv_emp    = FLD_C( 'conservative'        ,    'no'    ,     ''      ,         ''          ,   ''    )
213      sn_rcv_rnf    = FLD_C( 'coupled'             ,    'no'    ,     ''      ,         ''          ,   ''    )
214      sn_rcv_cal    = FLD_C( 'coupled'             ,    'no'    ,     ''      ,         ''          ,   ''    )
215      sn_rcv_iceflx = FLD_C( 'none'                ,    'no'    ,     ''      ,         ''          ,   ''    )
216      sn_rcv_co2    = FLD_C( 'none'                ,    'no'    ,     ''      ,         ''          ,   ''    )
217
218      REWIND( numnam )                    ! ... read namlist namsbc_cpl
219      READ  ( numnam, namsbc_cpl )
220
221      IF(lwp) THEN                        ! control print
222         WRITE(numout,*)
223         WRITE(numout,*)'sbc_cpl_init : namsbc_cpl namelist '
224         WRITE(numout,*)'~~~~~~~~~~~~'
225         WRITE(numout,*)'  received fields (mutiple ice categogies)'
226         WRITE(numout,*)'      10m wind module                 = ', TRIM(sn_rcv_w10m%cldes  ), ' (', TRIM(sn_rcv_w10m%clcat  ), ')'
227         WRITE(numout,*)'      stress module                   = ', TRIM(sn_rcv_taumod%cldes), ' (', TRIM(sn_rcv_taumod%clcat), ')'
228         WRITE(numout,*)'      surface stress                  = ', TRIM(sn_rcv_tau%cldes   ), ' (', TRIM(sn_rcv_tau%clcat   ), ')'
229         WRITE(numout,*)'                     - referential    = ', sn_rcv_tau%clvref
230         WRITE(numout,*)'                     - orientation    = ', sn_rcv_tau%clvor
231         WRITE(numout,*)'                     - mesh           = ', sn_rcv_tau%clvgrd
232         WRITE(numout,*)'      non-solar heat flux sensitivity = ', TRIM(sn_rcv_dqnsdt%cldes), ' (', TRIM(sn_rcv_dqnsdt%clcat), ')'
233         WRITE(numout,*)'      solar heat flux                 = ', TRIM(sn_rcv_qsr%cldes   ), ' (', TRIM(sn_rcv_qsr%clcat   ), ')'
234         WRITE(numout,*)'      non-solar heat flux             = ', TRIM(sn_rcv_qns%cldes   ), ' (', TRIM(sn_rcv_qns%clcat   ), ')'
235         WRITE(numout,*)'      freshwater budget               = ', TRIM(sn_rcv_emp%cldes   ), ' (', TRIM(sn_rcv_emp%clcat   ), ')'
236         WRITE(numout,*)'      runoffs                         = ', TRIM(sn_rcv_rnf%cldes   ), ' (', TRIM(sn_rcv_rnf%clcat   ), ')'
237         WRITE(numout,*)'      calving                         = ', TRIM(sn_rcv_cal%cldes   ), ' (', TRIM(sn_rcv_cal%clcat   ), ')'
238         WRITE(numout,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')'
239         WRITE(numout,*)'      atm co2                         = ', TRIM(sn_rcv_co2%cldes   ), ' (', TRIM(sn_rcv_co2%clcat   ), ')'
240         WRITE(numout,*)'  sent fields (multiple ice categories)'
241         WRITE(numout,*)'      surface temperature             = ', TRIM(sn_snd_temp%cldes  ), ' (', TRIM(sn_snd_temp%clcat  ), ')'
242         WRITE(numout,*)'      albedo                          = ', TRIM(sn_snd_alb%cldes   ), ' (', TRIM(sn_snd_alb%clcat   ), ')'
243         WRITE(numout,*)'      ice/snow thickness              = ', TRIM(sn_snd_thick%cldes ), ' (', TRIM(sn_snd_thick%clcat ), ')'
244         WRITE(numout,*)'      surface current                 = ', TRIM(sn_snd_crt%cldes   ), ' (', TRIM(sn_snd_crt%clcat   ), ')'
245         WRITE(numout,*)'                      - referential   = ', sn_snd_crt%clvref 
246         WRITE(numout,*)'                      - orientation   = ', sn_snd_crt%clvor
247         WRITE(numout,*)'                      - mesh          = ', sn_snd_crt%clvgrd
248         WRITE(numout,*)'      oce co2 flux                    = ', TRIM(sn_snd_co2%cldes   ), ' (', TRIM(sn_snd_co2%clcat   ), ')'
249      ENDIF
250
251      !                                   ! allocate sbccpl arrays
252      IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )
253     
254      ! ================================ !
255      !   Define the receive interface   !
256      ! ================================ !
257
258      ! Vectors: change of sign at north fold ONLY if on the local grid
259      IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' )   sf(jp_otx1:jp_itz2)%nsgn = -1.
260     
261      !
262      IF( TRIM( sn_rcv_tau%clvor  ) == 'local grid' ) THEN        ! already on local grid -> no need of the second grid
263            sf(jp_oty1)%clvgrd = sf(jp_oty2)%clvgrd   ! not needed but cleaner...
264            sf(jp_ity1)%clvgrd = sf(jp_ity2)%clvgrd   ! not needed but cleaner...
265      ENDIF
266      !
267
268      ! ================================ !
269      !     Define the send interface    !
270      ! ================================ !
271      ! for each field: define the OASIS name                           (ssnd(:)%clname)
272      !                 define send or not from the namelist parameters (ssnd(:)%laction)
273      !                 define the north fold type of lbc               (ssnd(:)%nsgn)
274     
275      ! default definitions of nsnd
276      ssnd(:)%laction = .FALSE.   ;   ssnd(:)%clgrid = 'T'   ;   ssnd(:)%nsgn = 1.  ; ssnd(:)%nct = 1
277         
278      !                                                      ! ------------------------- !
279      !                                                      !    Surface temperature    !
280      !                                                      ! ------------------------- !
281      ssnd(jps_toce)%clname = 'O_SSTSST'
282      ssnd(jps_tice)%clname = 'O_TepIce'
283      ssnd(jps_tmix)%clname = 'O_TepMix'
284      SELECT CASE( TRIM( sn_snd_temp%cldes ) )
285      CASE( 'none'         )       ! nothing to do
286      CASE( 'oce only'             )   ;   ssnd(   jps_toce             )%laction = .TRUE.
287      CASE( 'weighted oce and ice' )
288         ssnd( (/jps_toce, jps_tice/) )%laction = .TRUE.
289         IF ( TRIM( sn_snd_temp%clcat ) == 'yes' )  ssnd(jps_tice)%nct = jpl
290      CASE( 'mixed oce-ice'        )   ;   ssnd(   jps_tmix             )%laction = .TRUE.
291      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_temp%cldes' )
292      END SELECT
293     
294      !                                                      ! ------------------------- !
295      !                                                      !          Albedo           !
296      !                                                      ! ------------------------- !
297      ssnd(jps_albice)%clname = 'O_AlbIce' 
298      ssnd(jps_albmix)%clname = 'O_AlbMix'
299      SELECT CASE( TRIM( sn_snd_alb%cldes ) )
300      CASE( 'none'          )       ! nothing to do
301      CASE( 'weighted ice'  )   ;   ssnd(jps_albice)%laction = .TRUE.
302      CASE( 'mixed oce-ice' )   ;   ssnd(jps_albmix)%laction = .TRUE.
303      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_alb%cldes' )
304      END SELECT
305      !
306      ! Need to calculate oceanic albedo if
307      !     1. sending mixed oce-ice albedo or
308      !     2. receiving mixed oce-ice solar radiation
309      IF ( TRIM ( sn_snd_alb%cldes ) == 'mixed oce-ice' .OR. TRIM ( sn_rcv_qsr%cldes ) == 'mixed oce-ice' ) THEN
310         CALL albedo_oce( zaos, zacs )
311         ! Due to lack of information on nebulosity : mean clear/overcast sky
312         albedo_oce_mix(:,:) = ( zacs(:,:) + zaos(:,:) ) * 0.5
313      ENDIF
314
315      !                                                      ! ------------------------- !
316      !                                                      !  Ice fraction & Thickness !
317      !                                                      ! ------------------------- !
318      ssnd(jps_fice)%clname = 'OIceFrc'
319      ssnd(jps_hice)%clname = 'OIceTck'
320      ssnd(jps_hsnw)%clname = 'OSnwTck'
321      IF( k_ice /= 0 ) THEN
322         ssnd(jps_fice)%laction = .TRUE.                  ! if ice treated in the ocean (even in climato case)
323! Currently no namelist entry to determine sending of multi-category ice fraction so use the thickness entry for now
324         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_fice)%nct = jpl
325      ENDIF
326
327      SELECT CASE ( TRIM( sn_snd_thick%cldes ) )
328      CASE( 'none'         )       ! nothing to do
329      CASE( 'ice and snow' ) 
330         ssnd(jps_hice:jps_hsnw)%laction = .TRUE.
331         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) THEN
332            ssnd(jps_hice:jps_hsnw)%nct = jpl
333         ELSE
334            IF ( jpl > 1 ) THEN
335CALL ctl_stop( 'sbc_cpl_init: use weighted ice and snow option for sn_snd_thick%cldes if not exchanging category fields' )
336            ENDIF
337         ENDIF
338      CASE ( 'weighted ice and snow' ) 
339         ssnd(jps_hice:jps_hsnw)%laction = .TRUE.
340         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_hice:jps_hsnw)%nct = jpl
341      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_thick%cldes' )
342      END SELECT
343
344      !                                                      ! ------------------------- !
345      !                                                      !      Surface current      !
346      !                                                      ! ------------------------- !
347      !        ocean currents              !            ice velocities
348      ssnd(jps_ocx1)%clname = 'O_OCurx1'   ;   ssnd(jps_ivx1)%clname = 'O_IVelx1'
349      ssnd(jps_ocy1)%clname = 'O_OCury1'   ;   ssnd(jps_ivy1)%clname = 'O_IVely1'
350      ssnd(jps_ocz1)%clname = 'O_OCurz1'   ;   ssnd(jps_ivz1)%clname = 'O_IVelz1'
351      !
352      ssnd(jps_ocx1:jps_ivz1)%nsgn = -1.   ! vectors: change of the sign at the north fold
353
354      IF( sn_snd_crt%clvgrd == 'U,V' ) THEN
355         ssnd(jps_ocx1)%clgrid = 'U' ; ssnd(jps_ocy1)%clgrid = 'V'
356      ELSE IF( sn_snd_crt%clvgrd /= 'T' ) THEN 
357         CALL ctl_stop( 'sn_snd_crt%clvgrd must be equal to T' )
358         ssnd(jps_ocx1:jps_ivz1)%clgrid  = 'T'      ! all oce and ice components on the same unique grid
359      ENDIF
360      ssnd(jps_ocx1:jps_ivz1)%laction = .TRUE.   ! default: all are send
361      IF( TRIM( sn_snd_crt%clvref ) == 'spherical' )   ssnd( (/jps_ocz1, jps_ivz1/) )%laction = .FALSE. 
362      IF( TRIM( sn_snd_crt%clvor ) == 'eastward-northward' ) ssnd(jps_ocx1:jps_ivz1)%nsgn = 1.
363      SELECT CASE( TRIM( sn_snd_crt%cldes ) )
364      CASE( 'none'                 )   ;   ssnd(jps_ocx1:jps_ivz1)%laction = .FALSE.
365      CASE( 'oce only'             )   ;   ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE.
366      CASE( 'weighted oce and ice' )   !   nothing to do
367      CASE( 'mixed oce-ice'        )   ;   ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE.
368      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_crt%cldes' )
369      END SELECT
370
371      !                                                      ! ------------------------- !
372      !                                                      !          CO2 flux         !
373      !                                                      ! ------------------------- !
374      ssnd(jps_co2)%clname = 'O_CO2FLX' ;  IF( TRIM(sn_snd_co2%cldes) == 'coupled' )    ssnd(jps_co2 )%laction = .TRUE.
375      !
376      ! ================================ !
377      !   initialisation of the coupler  !
378      ! ================================ !
379
380      CALL cpl_prism_define(jpfld, jpsnd, sf)           
381      !
382      IF( ln_dm2dc .AND. ( cpl_prism_freq( jp_qsroce ) + cpl_prism_freq( jp_qsrmix ) /= 86400 ) )   &
383         &   CALL ctl_stop( 'sbc_cpl_init: diurnal cycle reconstruction (ln_dm2dc) needs daily couping for solar radiation' )
384
385      CALL wrk_dealloc( jpi,jpj, zacs, zaos )
386      !
387      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_init')
388      !
389   END SUBROUTINE sbc_cpl_init
390
391
392   SUBROUTINE sbc_cpl_rcv( kt, k_fsbc, k_ice )     
393      !!----------------------------------------------------------------------
394      !!             ***  ROUTINE sbc_cpl_rcv  ***
395      !!
396      !! ** Purpose :   provide the stress over the ocean and, if no sea-ice,
397      !!                provide the ocean heat and freshwater fluxes.
398      !!
399      !! ** Method  : - Receive all the atmospheric fields (stored in sf array). called at each time step.
400      !!                OASIS controls if there is something do receive or not. ninfo contains the info
401      !!                to know if the field was really received or not
402      !!
403      !!              --> If ocean stress was really received:
404      !!
405      !!                  - transform the received ocean stress vector from the received
406      !!                 referential and grid into an atmosphere-ocean stress in
407      !!                 the (i,j) ocean referencial and at the ocean velocity point.
408      !!                    The received stress are :
409      !!                     - defined by 3 components (if cartesian coordinate)
410      !!                            or by 2 components (if spherical)
411      !!                     - oriented along geographical   coordinate (if eastward-northward)
412      !!                            or  along the local grid coordinate (if local grid)
413      !!                     - given at U- and V-point, resp.   if received on 2 grids
414      !!                            or at T-point               if received on 1 grid
415      !!                    Therefore and if necessary, they are successively
416      !!                  processed in order to obtain them
417      !!                     first  as  2 components on the sphere
418      !!                     second as  2 components oriented along the local grid
419      !!                     third  as  2 components on the U,V grid
420      !!
421      !!              -->
422      !!
423      !!              - In 'ocean only' case, non solar and solar ocean heat fluxes
424      !!             and total ocean freshwater fluxes 
425      !!
426      !! ** Method  :   receive all fields from the atmosphere and transform
427      !!              them into ocean surface boundary condition fields
428      !!
429      !! ** Action  :   update  utau, vtau   ocean stress at U,V grid
430      !!                        taum, wndm   wind stres and wind speed module at T-point
431      !!                        qns          non solar heat fluxes including emp heat content    (ocean only case)
432      !!                                     and the latent heat flux of solid precip. melting
433      !!                        qsr          solar ocean heat fluxes   (ocean only case)
434      !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case)
435      !!----------------------------------------------------------------------
436      INTEGER, INTENT(in) ::   kt       ! ocean model time step index
437      INTEGER, INTENT(in) ::   k_fsbc   ! frequency of sbc (-> ice model) computation
438      INTEGER, INTENT(in) ::   k_ice    ! ice management in the sbc (=0/1/2/3)
439      !!
440      LOGICAL ::    llnewtx, llnewtau      ! update wind stress components and module??
441      INTEGER  ::   ji, jj, jn             ! dummy loop indices
442      REAL(wp) ::   zcumulneg, zcumulpos   ! temporary scalars     
443      REAL(wp) ::   zcoef                  ! temporary scalar
444      REAL(wp) ::   zrhoa  = 1.22          ! Air density kg/m3
445      REAL(wp) ::   zcdrag = 1.5e-3        ! drag coefficient
446      REAL(wp) ::   zzx, zzy               ! temporary variables
447      REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty 
448      !!----------------------------------------------------------------------
449      !
450      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_rcv')
451      !
452      CALL wrk_alloc( jpi,jpj, ztx, zty )
453
454      IF( kt == nit000 )   CALL sbc_cpl_init( k_ice )          ! initialisation
455
456      !                                                 ! Receive all the atmos. fields (including ice information)
457      !                                                      ! ========================= !
458      IF( sf(jp_otx1)%loasis ) THEN                      !  ocean stress components  !
459         !                                                   ! ========================= !
460         ! define sf(jp_otx1)%fnow(:,:,1) and sf(jp_oty1)%fnow(:,:,1): stress at U/V point along model grid
461         ! => need to be done only when we receive the field
462         IF(  sf(jp_otx1)%ninfo == OASIS_Rcv ) THEN
463            !
464            IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN            ! 2 components on the sphere
465               !                                                       ! (cartesian to spherical -> 3 to 2 components)
466               !
467               CALL geo2oce( sf(jp_otx1)%fnow(:,:,1), sf(jp_oty1)%fnow(:,:,1), sf(jp_otz1)%fnow(:,:,1),   &
468                  &          sf(jp_otx1)%clvgrd, ztx, zty )
469               sf(jp_otx1)%fnow(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 1st grid
470               sf(jp_oty1)%fnow(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 1st grid
471               !
472               IF( sf(jp_otx2)%loasis ) THEN
473                  CALL geo2oce( sf(jp_otx2)%fnow(:,:,1), sf(jp_oty2)%fnow(:,:,1), sf(jp_otz2)%fnow(:,:,1),   &
474                     &          sf(jp_otx2)%clvgrd, ztx, zty )
475                  sf(jp_otx2)%fnow(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 2nd grid
476                  sf(jp_oty2)%fnow(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 2nd grid
477               ENDIF
478               !
479            ENDIF
480            !
481            IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN   ! 2 components oriented along the local grid
482               !                                                       ! (geographical to local grid -> rotate the components)
483               CALL rot_rep( sf(jp_otx1)%fnow(:,:,1), sf(jp_oty1)%fnow(:,:,1), sf(jp_otx1)%clvgrd, 'en->i', ztx )   
484               IF( sf(jp_otx2)%loasis ) THEN
485                  CALL rot_rep( sf(jp_otx2)%fnow(:,:,1), sf(jp_oty2)%fnow(:,:,1), sf(jp_otx2)%clvgrd, 'en->j', zty )   
486               ELSE 
487                  CALL rot_rep( sf(jp_otx1)%fnow(:,:,1), sf(jp_oty1)%fnow(:,:,1), sf(jp_otx1)%clvgrd, 'en->j', zty ) 
488               ENDIF
489               sf(jp_otx1)%fnow(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid
490               sf(jp_oty1)%fnow(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 2nd grid
491            ENDIF
492            !                             
493            IF( sf(jp_otx1)%clvgrd == 'T' ) THEN
494               DO jj = 2, jpjm1                                          ! T ==> (U,V)
495                  DO ji = fs_2, fs_jpim1   ! vector opt.
496                     sf(jp_otx1)%fnow(ji,jj,1) = 0.5 * ( sf(jp_otx1)%fnow(ji+1,jj  ,1) + sf(jp_otx1)%fnow(ji,jj,1) )
497                     sf(jp_oty1)%fnow(ji,jj,1) = 0.5 * ( sf(jp_oty1)%fnow(ji  ,jj+1,1) + sf(jp_oty1)%fnow(ji,jj,1) )
498                  END DO
499               END DO
500               CALL lbc_lnk( sf(jp_otx1)%fnow(:,:,1), 'U',  -1. )   ;   CALL lbc_lnk( sf(jp_oty1)%fnow(:,:,1), 'V',  -1. )
501            ENDIF
502            llnewtx = .TRUE.
503         ELSE
504            llnewtx = .FALSE.
505         ENDIF
506         !                                                   ! ========================= !
507      ELSE                                                   !   No dynamical coupling   !
508         !                                                   ! ========================= !
509         sf(jp_otx1)%fnow(:,:,1) = 0.e0                               ! here simply set to zero
510         sf(jp_oty1)%fnow(:,:,1) = 0.e0                               ! an external read in a file can be added instead
511         llnewtx = .TRUE.
512         !
513      ENDIF
514     
515      !                                                      ! ========================= !
516      !                                                      !    wind stress module     !   (taum)
517      !                                                      ! ========================= !
518      !
519      IF( .NOT. sf(jp_taum)%loasis ) THEN                    ! compute wind stress module from its components if not received
520         ! => need to be done only when otx1 was changed
521         IF( llnewtx ) THEN
522!CDIR NOVERRCHK
523            DO jj = 2, jpjm1
524!CDIR NOVERRCHK
525               DO ji = fs_2, fs_jpim1   ! vect. opt.
526                  zzx = sf(jp_otx1)%fnow(ji-1,jj  ,1) + sf(jp_otx1)%fnow(ji,jj,1)
527                  zzy = sf(jp_oty1)%fnow(ji  ,jj-1,1) + sf(jp_oty1)%fnow(ji,jj,1)
528                  sf(jp_taum)%fnow(ji,jj,1) = 0.5 * SQRT( zzx * zzx + zzy * zzy )
529               END DO
530            END DO
531            CALL lbc_lnk( sf(jp_taum)%fnow(:,:,1), 'T', 1. )
532            llnewtau = .TRUE.
533         ELSE
534            llnewtau = .FALSE.
535         ENDIF
536      ELSE
537         llnewtau = sf(jp_taum)%ninfo == OASIS_Rcv
538         ! Stress module can be negative when received (interpolation problem)
539         IF( llnewtau ) THEN
540            sf(jp_taum)%fnow(:,:,1) = MAX( 0._wp, sf(jp_taum)%fnow(:,:,1) )
541         ENDIF
542      ENDIF
543     
544      !                                                      ! ========================= !
545      !                                                      !      10 m wind speed      !   (wndm)
546      !                                                      ! ========================= !
547      !
548      IF( .NOT. sf(jp_w10m)%loasis ) THEN                    ! compute wind spreed from wind stress module if not received 
549         ! => need to be done only when taumod was changed
550         IF( llnewtau ) THEN
551            zcoef = 1. / ( zrhoa * zcdrag ) 
552!CDIR NOVERRCHK
553            DO jj = 1, jpj
554!CDIR NOVERRCHK
555               DO ji = 1, jpi 
556                  wndm(ji,jj) = SQRT( sf(jp_taum)%fnow(ji,jj,1) * zcoef )
557               END DO
558            END DO
559         ENDIF
560      ELSE
561         IF ( sf(jp_w10m)%ninfo == OASIS_Rcv ) wndm(:,:) = sf(jp_w10m)%fnow(:,:,1)
562      ENDIF
563
564      ! u(v)tau and taum will be modified by ice model
565      ! -> need to be reset before each call of the ice/fsbc     
566      IF( MOD( kt-1, k_fsbc ) == 0 ) THEN
567         !
568         utau(:,:) = sf(jp_otx1)%fnow(:,:,1)
569         vtau(:,:) = sf(jp_oty1)%fnow(:,:,1)
570         taum(:,:) = sf(jp_taum)%fnow(:,:,1)
571         CALL iom_put( "taum_oce", taum )   ! output wind stress module
572         
573      ENDIF
574
575#if defined key_cpl_carbon_cycle
576      !                                                              ! atmosph. CO2 (ppm)
577      IF( sf(jp_co2)%loasis )   atm_co2(:,:) = sf(jp_co2)%fnow(:,:,1)
578#endif
579
580      !                                                      ! ========================= !
581      IF( k_ice <= 1 ) THEN                                  !  heat & freshwater fluxes ! (Ocean only case)
582         !                                                   ! ========================= !
583         !
584         !                                                       ! total freshwater fluxes over the ocean (emp)
585         SELECT CASE( TRIM( sn_rcv_emp%cldes ) )                                    ! evaporation - precipitation
586         CASE( 'conservative' )
587            emp(:,:) = sf(jp_tevp)%fnow(:,:,1) - ( sf(jp_rain)%fnow(:,:,1) + sf(jp_snow)%fnow(:,:,1) )
588         CASE( 'oce only', 'oce and ice' )
589            emp(:,:) = sf(jp_oemp)%fnow(:,:,1)
590         CASE default
591            CALL ctl_stop( 'sbc_cpl_rcv: wrong definition of sn_rcv_emp%cldes' )
592         END SELECT
593         !
594         !                                                        ! runoffs and calving (added in emp)
595         IF( sf(jp_rnf)%loasis )   emp(:,:) = emp(:,:) - sf(jp_rnf)%fnow(:,:,1)
596         IF( sf(jp_cal)%loasis )   emp(:,:) = emp(:,:) - sf(jp_cal)%fnow(:,:,1)
597         !
598!!gm :  this seems to be internal cooking, not sure to need that in a generic interface
599!!gm                                       at least should be optional...
600!!         IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' ) THEN     ! add to the total freshwater budget
601!!            ! remove negative runoff
602!!            zcumulpos = SUM( MAX( sf(jp_rnf)%fnow(:,:,1), 0.e0 ) * e1t(:,:) * e2t(:,:) * tmask_i(:,:) )
603!!            zcumulneg = SUM( MIN( sf(jp_rnf)%fnow(:,:,1), 0.e0 ) * e1t(:,:) * e2t(:,:) * tmask_i(:,:) )
604!!            IF( lk_mpp )   CALL mpp_sum( zcumulpos )   ! sum over the global domain
605!!            IF( lk_mpp )   CALL mpp_sum( zcumulneg )
606!!            IF( zcumulpos /= 0. ) THEN                 ! distribute negative runoff on positive runoff grid points
607!!               zcumulneg = 1.e0 + zcumulneg / zcumulpos
608!!               sf(jp_rnf)%fnow(:,:,1) = MAX( sf(jp_rnf)%fnow(:,:,1), 0.e0 ) * zcumulneg
609!!            ENDIF     
610!!            ! add runoff to e-p
611!!            emp(:,:) = emp(:,:) - sf(jp_rnf)%fnow(:,:,1)
612!!         ENDIF
613!!gm  end of internal cooking
614         !
615         !                                                       ! non solar heat flux over the ocean (qns)
616         IF( sf(jp_qnsoce)%loasis )   qns(:,:) = sf(jp_qnsoce)%fnow(:,:,1)
617         IF( sf(jp_qnsmix)%loasis )   qns(:,:) = sf(jp_qnsmix)%fnow(:,:,1)
618         ! add the latent heat of solid precip. melting
619         IF( sf(jp_snow  )%loasis )   THEN                         ! update qns over the free ocean with:
620              qns(:,:) = qns(:,:) - sf(jp_snow)%fnow(:,:,1) * lfus  & ! energy for melting solid precipitation over the free ocean
621           &           - emp(:,:) * sst_m(:,:) * rcp                   ! remove heat content due to mass flux (assumed to be at SST)
622         ENDIF
623
624         !                                                       ! solar flux over the ocean          (qsr)
625         IF( sf(jp_qsroce)%loasis )   qsr(:,:) = sf(jp_qsroce)%fnow(:,:,1)
626         IF( sf(jp_qsrmix)%loasis )   qsr(:,:) = sf(jp_qsrmix)%fnow(:,:,1)
627         IF( ln_dm2dc )   qsr(:,:) = sbc_dcy( qsr )                           ! modify qsr to include the diurnal cycle
628         !
629 
630      ENDIF
631      !
632      CALL wrk_dealloc( jpi,jpj, ztx, zty )
633      !
634      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_rcv')
635      !
636   END SUBROUTINE sbc_cpl_rcv
637   
638
639   SUBROUTINE sbc_cpl_ice_tau( p_taui, p_tauj )     
640      !!----------------------------------------------------------------------
641      !!             ***  ROUTINE sbc_cpl_ice_tau  ***
642      !!
643      !! ** Purpose :   provide the stress over sea-ice in coupled mode
644      !!
645      !! ** Method  :   transform the received stress from the atmosphere into
646      !!             an atmosphere-ice stress in the (i,j) ocean referencial
647      !!             and at the velocity point of the sea-ice model (cp_ice_msh):
648      !!                'C'-grid : i- (j-) components given at U- (V-) point
649      !!                'I'-grid : B-grid lower-left corner: both components given at I-point
650      !!
651      !!                The received stress are :
652      !!                 - defined by 3 components (if cartesian coordinate)
653      !!                        or by 2 components (if spherical)
654      !!                 - oriented along geographical   coordinate (if eastward-northward)
655      !!                        or  along the local grid coordinate (if local grid)
656      !!                 - given at U- and V-point, resp.   if received on 2 grids
657      !!                        or at a same point (T or I) if received on 1 grid
658      !!                Therefore and if necessary, they are successively
659      !!             processed in order to obtain them
660      !!                 first  as  2 components on the sphere
661      !!                 second as  2 components oriented along the local grid
662      !!                 third  as  2 components on the cp_ice_msh point
663      !!
664      !!                In 'oce and ice' case, only one vector stress field
665      !!             is received. It has already been processed in sbc_cpl_rcv
666      !!             so that it is now defined as (i,j) components given at U-
667      !!             and V-points, respectively. Therefore, here only the third
668      !!             transformation is done and only if the ice-grid is a 'I'-grid.
669      !!
670      !! ** Action  :   return ptau_i, ptau_j, the stress over the ice at cp_ice_msh point
671      !!----------------------------------------------------------------------
672      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_taui   ! i- & j-components of atmos-ice stress [N/m2]
673      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_tauj   ! at I-point (B-grid) or U & V-point (C-grid)
674      !!
675      INTEGER ::   ji, jj                          ! dummy loop indices
676      INTEGER ::   itx                             ! index of taux over ice
677      REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty 
678      !!----------------------------------------------------------------------
679      !
680      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_tau')
681      !
682      CALL wrk_alloc( jpi,jpj, ztx, zty )
683
684      IF( sf(jp_itx1)%loasis ) THEN   ;   itx =  jp_itx1   
685      ELSE                                ;   itx =  jp_otx1
686      ENDIF
687
688      ! do something only if we just received the stress from atmosphere
689      IF(  sf(itx)%ninfo == OASIS_Rcv ) THEN
690
691         !                                                      ! ======================= !
692         IF( sf(jp_itx1)%loasis ) THEN                      !   ice stress received   !
693            !                                                   ! ======================= !
694           
695            IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN            ! 2 components on the sphere
696               !                                                       ! (cartesian to spherical -> 3 to 2 components)
697               CALL geo2oce(  sf(jp_itx1)%fnow(:,:,1), sf(jp_ity1)%fnow(:,:,1), sf(jp_itz1)%fnow(:,:,1),   &
698                  &          sf(jp_itx1)%clvgrd, ztx, zty )
699               sf(jp_itx1)%fnow(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 1st grid
700               sf(jp_ity1)%fnow(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 1st grid
701               !
702               IF( sf(jp_itx2)%loasis ) THEN
703                  CALL geo2oce( sf(jp_itx2)%fnow(:,:,1), sf(jp_ity2)%fnow(:,:,1), sf(jp_itz2)%fnow(:,:,1),   &
704                     &          sf(jp_itx2)%clvgrd, ztx, zty )
705                  sf(jp_itx2)%fnow(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 2nd grid
706                  sf(jp_ity2)%fnow(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 2nd grid
707               ENDIF
708               !
709            ENDIF
710            !
711            IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN   ! 2 components oriented along the local grid
712               !                                                       ! (geographical to local grid -> rotate the components)
713               CALL rot_rep( sf(jp_itx1)%fnow(:,:,1), sf(jp_ity1)%fnow(:,:,1), sf(jp_itx1)%clvgrd, 'en->i', ztx )   
714               IF( sf(jp_itx2)%loasis ) THEN
715                  CALL rot_rep( sf(jp_itx2)%fnow(:,:,1), sf(jp_ity2)%fnow(:,:,1), sf(jp_itx2)%clvgrd, 'en->j', zty )   
716               ELSE
717                  CALL rot_rep( sf(jp_itx1)%fnow(:,:,1), sf(jp_ity1)%fnow(:,:,1), sf(jp_itx1)%clvgrd, 'en->j', zty ) 
718               ENDIF
719               sf(jp_itx1)%fnow(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid
720               sf(jp_ity1)%fnow(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 1st grid
721            ENDIF
722            !                                                   ! ======================= !
723         ELSE                                                   !     use ocean stress    !
724            !                                                   ! ======================= !
725            sf(jp_itx1)%fnow(:,:,1) = sf(jp_otx1)%fnow(:,:,1)
726            sf(jp_ity1)%fnow(:,:,1) = sf(jp_oty1)%fnow(:,:,1)
727            !
728         ENDIF
729
730         !                                                      ! ======================= !
731         !                                                      !     put on ice grid     !
732         !                                                      ! ======================= !
733         !   
734         !                                                  j+1   j     -----V---F
735         ! ice stress on ice velocity point (cp_ice_msh)                 !       |
736         ! (C-grid ==>(U,V) or B-grid ==> I or F)                 j      |   T   U
737         !                                                               |       |
738         !                                                   j    j-1   -I-------|
739         !                                               (for I)         |       |
740         !                                                              i-1  i   i
741         !                                                               i      i+1 (for I)
742         SELECT CASE ( cp_ice_msh )
743            !
744         CASE( 'I' )                                         ! B-grid ==> I
745            SELECT CASE ( sf(jp_itx1)%clvgrd )
746            CASE( 'U' )
747               DO jj = 2, jpjm1                                   ! (U,V) ==> I
748                  DO ji = 2, jpim1   ! NO vector opt.
749                     p_taui(ji,jj) = 0.5 * ( sf(jp_itx1)%fnow(ji-1,jj  ,1) + sf(jp_itx1)%fnow(ji-1,jj-1,1) )
750                     p_tauj(ji,jj) = 0.5 * ( sf(jp_ity1)%fnow(ji  ,jj-1,1) + sf(jp_ity1)%fnow(ji-1,jj-1,1) )
751                  END DO
752               END DO
753            CASE( 'F' )
754               DO jj = 2, jpjm1                                   ! F ==> I
755                  DO ji = 2, jpim1   ! NO vector opt.
756                     p_taui(ji,jj) = sf(jp_itx1)%fnow(ji-1,jj-1,1)
757                     p_tauj(ji,jj) = sf(jp_ity1)%fnow(ji-1,jj-1,1)
758                  END DO
759               END DO
760            CASE( 'T' )
761               DO jj = 2, jpjm1                                   ! T ==> I
762                  DO ji = 2, jpim1   ! NO vector opt.
763                     p_taui(ji,jj) = 0.25 * ( sf(jp_itx1)%fnow(ji,jj  ,1) + sf(jp_itx1)%fnow(ji-1,jj  ,1)   &
764                        &                   + sf(jp_itx1)%fnow(ji,jj-1,1) + sf(jp_itx1)%fnow(ji-1,jj-1,1) ) 
765                     p_tauj(ji,jj) = 0.25 * ( sf(jp_ity1)%fnow(ji,jj  ,1) + sf(jp_ity1)%fnow(ji-1,jj  ,1)   &
766                        &                   + sf(jp_oty1)%fnow(ji,jj-1,1) + sf(jp_ity1)%fnow(ji-1,jj-1,1) )
767                  END DO
768               END DO
769            CASE( 'I' )
770               p_taui(:,:) = sf(jp_itx1)%fnow(:,:,1)                   ! I ==> I
771               p_tauj(:,:) = sf(jp_ity1)%fnow(:,:,1)
772            END SELECT
773            IF( sf(jp_itx1)%clvgrd /= 'I' ) THEN
774               CALL lbc_lnk( p_taui, 'I',  -1. )   ;   CALL lbc_lnk( p_tauj, 'I',  -1. )
775            ENDIF
776            !
777         CASE( 'F' )                                         ! B-grid ==> F
778            SELECT CASE ( sf(jp_itx1)%clvgrd )
779            CASE( 'U' )
780               DO jj = 2, jpjm1                                   ! (U,V) ==> F
781                  DO ji = fs_2, fs_jpim1   ! vector opt.
782                     p_taui(ji,jj) = 0.5 * ( sf(jp_itx1)%fnow(ji,jj,1) + sf(jp_itx1)%fnow(ji  ,jj+1,1) )
783                     p_tauj(ji,jj) = 0.5 * ( sf(jp_ity1)%fnow(ji,jj,1) + sf(jp_ity1)%fnow(ji+1,jj  ,1) )
784                  END DO
785               END DO
786            CASE( 'I' )
787               DO jj = 2, jpjm1                                   ! I ==> F
788                  DO ji = 2, jpim1   ! NO vector opt.
789                     p_taui(ji,jj) = sf(jp_itx1)%fnow(ji+1,jj+1,1)
790                     p_tauj(ji,jj) = sf(jp_ity1)%fnow(ji+1,jj+1,1)
791                  END DO
792               END DO
793            CASE( 'T' )
794               DO jj = 2, jpjm1                                   ! T ==> F
795                  DO ji = 2, jpim1   ! NO vector opt.
796                     p_taui(ji,jj) = 0.25 * ( sf(jp_itx1)%fnow(ji,jj  ,1) + sf(jp_itx1)%fnow(ji+1,jj  ,1)   &
797                        &                   + sf(jp_itx1)%fnow(ji,jj+1,1) + sf(jp_itx1)%fnow(ji+1,jj+1,1) ) 
798                     p_tauj(ji,jj) = 0.25 * ( sf(jp_ity1)%fnow(ji,jj  ,1) + sf(jp_ity1)%fnow(ji+1,jj  ,1)   &
799                        &                   + sf(jp_ity1)%fnow(ji,jj+1,1) + sf(jp_ity1)%fnow(ji+1,jj+1,1) )
800                  END DO
801               END DO
802            CASE( 'F' )
803               p_taui(:,:) = sf(jp_itx1)%fnow(:,:,1)                   ! F ==> F
804               p_tauj(:,:) = sf(jp_ity1)%fnow(:,:,1)
805            END SELECT
806            IF( sf(jp_itx1)%clvgrd /= 'F' ) THEN
807               CALL lbc_lnk( p_taui, 'F',  -1. )   ;   CALL lbc_lnk( p_tauj, 'F',  -1. )
808            ENDIF
809            !
810         CASE( 'C' )                                         ! C-grid ==> U,V
811            SELECT CASE ( sf(jp_itx1)%clvgrd )
812            CASE( 'U' )
813               p_taui(:,:) = sf(jp_itx1)%fnow(:,:,1)                   ! (U,V) ==> (U,V)
814               p_tauj(:,:) = sf(jp_ity1)%fnow(:,:,1)
815            CASE( 'F' )
816               DO jj = 2, jpjm1                                   ! F ==> (U,V)
817                  DO ji = fs_2, fs_jpim1   ! vector opt.
818                     p_taui(ji,jj) = 0.5 * ( sf(jp_itx1)%fnow(ji,jj,1) + sf(jp_itx1)%fnow(ji  ,jj-1,1) )
819                     p_tauj(ji,jj) = 0.5 * ( sf(jp_ity1)%fnow(jj,jj,1) + sf(jp_ity1)%fnow(ji-1,jj  ,1) )
820                  END DO
821               END DO
822            CASE( 'T' )
823               DO jj = 2, jpjm1                                   ! T ==> (U,V)
824                  DO ji = fs_2, fs_jpim1   ! vector opt.
825                     p_taui(ji,jj) = 0.5 * ( sf(jp_itx1)%fnow(ji+1,jj  ,1) + sf(jp_itx1)%fnow(ji,jj,1) )
826                     p_tauj(ji,jj) = 0.5 * ( sf(jp_ity1)%fnow(ji  ,jj+1,1) + sf(jp_ity1)%fnow(ji,jj,1) )
827                  END DO
828               END DO
829            CASE( 'I' )
830               DO jj = 2, jpjm1                                   ! I ==> (U,V)
831                  DO ji = 2, jpim1   ! NO vector opt.
832                     p_taui(ji,jj) = 0.5 * ( sf(jp_itx1)%fnow(ji+1,jj+1,1) + sf(jp_itx1)%fnow(ji+1,jj  ,1) )
833                     p_tauj(ji,jj) = 0.5 * ( sf(jp_ity1)%fnow(ji+1,jj+1,1) + sf(jp_ity1)%fnow(ji  ,jj+1,1) )
834                  END DO
835               END DO
836            END SELECT
837            IF( sf(jp_itx1)%clvgrd /= 'U' ) THEN
838               CALL lbc_lnk( p_taui, 'U',  -1. )   ;   CALL lbc_lnk( p_tauj, 'V',  -1. )
839            ENDIF
840         END SELECT
841
842      ENDIF
843      !   
844      CALL wrk_dealloc( jpi,jpj, ztx, zty )
845      !
846      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_ice_tau')
847      !
848   END SUBROUTINE sbc_cpl_ice_tau
849   
850
851   SUBROUTINE sbc_cpl_ice_flx( p_frld  , palbi   , psst    , pist    )
852      !!----------------------------------------------------------------------
853      !!             ***  ROUTINE sbc_cpl_ice_flx  ***
854      !!
855      !! ** Purpose :   provide the heat and freshwater fluxes of the
856      !!              ocean-ice system.
857      !!
858      !! ** Method  :   transform the fields received from the atmosphere into
859      !!             surface heat and fresh water boundary condition for the
860      !!             ice-ocean system. The following fields are provided:
861      !!              * total non solar, solar and freshwater fluxes (qns_tot,
862      !!             qsr_tot and emp_tot) (total means weighted ice-ocean flux)
863      !!             NB: emp_tot include runoffs and calving.
864      !!              * fluxes over ice (qns_ice, qsr_ice, emp_ice) where
865      !!             emp_ice = sublimation - solid precipitation as liquid
866      !!             precipitation are re-routed directly to the ocean and
867      !!             runoffs and calving directly enter the ocean.
868      !!              * solid precipitation (sprecip), used to add to qns_tot
869      !!             the heat lost associated to melting solid precipitation
870      !!             over the ocean fraction.
871      !!       ===>> CAUTION here this changes the net heat flux received from
872      !!             the atmosphere
873      !!
874      !!                  - the fluxes have been separated from the stress as
875      !!                 (a) they are updated at each ice time step compare to
876      !!                 an update at each coupled time step for the stress, and
877      !!                 (b) the conservative computation of the fluxes over the
878      !!                 sea-ice area requires the knowledge of the ice fraction
879      !!                 after the ice advection and before the ice thermodynamics,
880      !!                 so that the stress is updated before the ice dynamics
881      !!                 while the fluxes are updated after it.
882      !!
883      !! ** Action  :   update at each nf_ice time step:
884      !!                   qns_tot, qsr_tot  non-solar and solar total heat fluxes
885      !!                   qns_ice, qsr_ice  non-solar and solar heat fluxes over the ice
886      !!                   emp_tot            total evaporation - precipitation(liquid and solid) (-runoff)(-calving)
887      !!                   emp_ice            ice sublimation - solid precipitation over the ice
888      !!                   dqns_ice           d(non-solar heat flux)/d(Temperature) over the ice
889      !!                   sprecip             solid precipitation over the ocean 
890      !!----------------------------------------------------------------------
891      REAL(wp), INTENT(in   ), DIMENSION(:,:)   ::   p_frld     ! lead fraction                [0 to 1]
892      ! optional arguments, used only in 'mixed oce-ice' case
893      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   palbi   ! ice albedo
894      REAL(wp), INTENT(in   ), DIMENSION(:,:  ), OPTIONAL ::   psst    ! sea surface temperature     [Celcius]
895      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   pist    ! ice surface temperature     [Kelvin]
896      !
897      INTEGER ::   jl   ! dummy loop index
898      REAL(wp), POINTER, DIMENSION(:,:) ::   zcptn, ztmp, zicefr
899      !!----------------------------------------------------------------------
900      !
901      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_flx')
902      !
903      CALL wrk_alloc( jpi,jpj, zcptn, ztmp, zicefr )
904
905      zicefr(:,:) = 1.- p_frld(:,:)
906      zcptn(:,:) = rcp * sst_m(:,:)
907      !
908      !                                                      ! ========================= !
909      !                                                      !    freshwater budget      !   (emp)
910      !                                                      ! ========================= !
911      !
912      !                                                           ! total Precipitations - total Evaporation (emp_tot)
913      !                                                           ! solid precipitation  - sublimation       (emp_ice)
914      !                                                           ! solid Precipitation                      (sprecip)
915      SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
916      CASE( 'conservative'  )   ! received fields: jp_rain, jp_snow, jp_ievp, jp_tevp
917         sprecip(:,:) = sf(jp_snow)%fnow(:,:,1)                 ! May need to ensure positive here
918         tprecip(:,:) = sf(jp_rain)%fnow(:,:,1) + sprecip (:,:) ! May need to ensure positive here
919         emp_tot(:,:) = sf(jp_tevp)%fnow(:,:,1) - tprecip(:,:)
920         emp_ice(:,:) = sf(jp_ievp)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1)
921                           CALL iom_put( 'rain'         , sf(jp_rain)%fnow(:,:,1)              )   ! liquid precipitation
922         IF( lk_diaar5 )   CALL iom_put( 'hflx_rain_cea', sf(jp_rain)%fnow(:,:,1) * zcptn(:,:) )   ! heat flux from liq. precip.
923         ztmp(:,:) = sf(jp_tevp)%fnow(:,:,1) - sf(jp_ievp)%fnow(:,:,1) * zicefr(:,:)
924                           CALL iom_put( 'evap_ao_cea'  , ztmp                            )   ! ice-free oce evap (cell average)
925         IF( lk_diaar5 )   CALL iom_put( 'hflx_evap_cea', ztmp(:,:         ) * zcptn(:,:) )   ! heat flux from from evap (cell ave)
926      CASE( 'oce and ice'   )   ! received fields: jp_sbpr, jp_semp, jp_oemp, jp_ievp
927         emp_tot(:,:) = p_frld(:,:) * sf(jp_oemp)%fnow(:,:,1) + zicefr(:,:) * sf(jp_sbpr)%fnow(:,:,1)
928         emp_ice(:,:) = sf(jp_semp)%fnow(:,:,1)
929         sprecip(:,:) = - sf(jp_semp)%fnow(:,:,1) + sf(jp_ievp)%fnow(:,:,1)
930      END SELECT
931
932      CALL iom_put( 'snowpre'    , sprecip                                )   ! Snow
933      CALL iom_put( 'snow_ao_cea', sprecip(:,:         ) * p_frld(:,:)    )   ! Snow        over ice-free ocean  (cell average)
934      CALL iom_put( 'snow_ai_cea', sprecip(:,:         ) * zicefr(:,:)    )   ! Snow        over sea-ice         (cell average)
935      CALL iom_put( 'subl_ai_cea', sf(jp_ievp)%fnow(:,:,1) * zicefr(:,:) )   ! Sublimation over sea-ice         (cell average)
936      !   
937      !                                                           ! runoffs and calving (put in emp_tot)
938      IF( sf(jp_rnf)%loasis ) THEN
939         emp_tot(:,:) = emp_tot(:,:) - sf(jp_rnf)%fnow(:,:,1)
940                           CALL iom_put( 'runoffs'      , sf(jp_rnf)%fnow(:,:,1)              )   ! rivers
941         IF( lk_diaar5 )   CALL iom_put( 'hflx_rnf_cea' , sf(jp_rnf)%fnow(:,:,1) * zcptn(:,:) )   ! heat flux from rivers
942      ENDIF
943      IF( sf(jp_cal)%loasis ) THEN
944         emp_tot(:,:) = emp_tot(:,:) - sf(jp_cal)%fnow(:,:,1)
945         CALL iom_put( 'calving', sf(jp_cal)%fnow(:,:,1) )
946      ENDIF
947      !
948!!gm :  this seems to be internal cooking, not sure to need that in a generic interface
949!!gm                                       at least should be optional...
950!!       ! remove negative runoff                            ! sum over the global domain
951!!       zcumulpos = SUM( MAX( sf(jp_rnf)%fnow(:,:,1), 0.e0 ) * e1t(:,:) * e2t(:,:) * tmask_i(:,:) )
952!!       zcumulneg = SUM( MIN( sf(jp_rnf)%fnow(:,:,1), 0.e0 ) * e1t(:,:) * e2t(:,:) * tmask_i(:,:) )
953!!       IF( lk_mpp )   CALL mpp_sum( zcumulpos )
954!!       IF( lk_mpp )   CALL mpp_sum( zcumulneg )
955!!       IF( zcumulpos /= 0. ) THEN                          ! distribute negative runoff on positive runoff grid points
956!!          zcumulneg = 1.e0 + zcumulneg / zcumulpos
957!!          sf(jp_rnf)%fnow(:,:,1) = MAX( sf(jp_rnf)%fnow(:,:,1), 0.e0 ) * zcumulneg
958!!       ENDIF     
959!!       emp_tot(:,:) = emp_tot(:,:) - sf(jp_rnf)%fnow(:,:,1)   ! add runoff to e-p
960!!
961!!gm  end of internal cooking
962
963      !                                                      ! ========================= !
964      SELECT CASE( TRIM( sn_rcv_qns%cldes ) )                !   non solar heat fluxes   !   (qns)
965      !                                                      ! ========================= !
966      CASE( 'oce only' )                                     ! the required field is directly provided
967         qns_tot(:,:  ) = sf(jp_qnsoce)%fnow(:,:,1)
968      CASE( 'conservative' )                                      ! the required fields are directly provided
969         qns_tot(:,:  ) = sf(jp_qnsmix)%fnow(:,:,1)
970         IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
971            qns_ice(:,:,1:jpl) = sf(jp_qnsice)%fnow(:,:,1:jpl)
972         ELSE
973            ! Set all category values equal for the moment
974            DO jl=1,jpl
975               qns_ice(:,:,jl) = sf(jp_qnsice)%fnow(:,:,1)
976            ENDDO
977         ENDIF
978      CASE( 'oce and ice' )       ! the total flux is computed from ocean and ice fluxes
979         qns_tot(:,:  ) =  p_frld(:,:) * sf(jp_qnsoce)%fnow(:,:,1)
980         IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
981            DO jl=1,jpl
982               qns_tot(:,:   ) = qns_tot(:,:) + a_i(:,:,jl) * sf(jp_qnsice)%fnow(:,:,jl)   
983               qns_ice(:,:,jl) = sf(jp_qnsice)%fnow(:,:,jl)
984            ENDDO
985         ELSE
986            DO jl=1,jpl
987               qns_tot(:,:   ) = qns_tot(:,:) + zicefr(:,:) * sf(jp_qnsice)%fnow(:,:,1)
988               qns_ice(:,:,jl) = sf(jp_qnsice)%fnow(:,:,1)
989            ENDDO
990         ENDIF
991      CASE( 'mixed oce-ice' )     ! the ice flux is cumputed from the total flux, the SST and ice informations
992! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
993         qns_tot(:,:  ) = sf(jp_qnsmix)%fnow(:,:,1)
994         qns_ice(:,:,1) = sf(jp_qnsmix)%fnow(:,:,1)    &
995            &            + sf(jp_dqnsdt)%fnow(:,:,1) * ( pist(:,:,1) - ( (rt0 + psst(:,:  ) ) * p_frld(:,:)   &
996            &                                                   +          pist(:,:,1)   * zicefr(:,:) ) )
997      END SELECT
998      ztmp(:,:) = p_frld(:,:) * sprecip(:,:) * lfus
999      qns_tot(:,:) = qns_tot(:,:)                         &            ! qns_tot update over free ocean with:
1000         &          - ztmp(:,:)                           &            ! remove the latent heat flux of solid precip. melting
1001         &          - (  emp_tot(:,:)                     &            ! remove the heat content of mass flux (assumed to be at SST)
1002         &             - emp_ice(:,:) * zicefr(:,:)  ) * zcptn(:,:) 
1003      IF( lk_diaar5 )   CALL iom_put( 'hflx_snow_cea', ztmp + sprecip(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average)
1004!!gm
1005!!    currently it is taken into account in leads budget but not in the qns_tot, and thus not in
1006!!    the flux that enter the ocean....
1007!!    moreover 1 - it is not diagnose anywhere....
1008!!             2 - it is unclear for me whether this heat lost is taken into account in the atmosphere or not...
1009!!
1010!! similar job should be done for snow and precipitation temperature
1011      !                                     
1012      IF( sf(jp_cal)%loasis ) THEN                            ! Iceberg melting
1013         ztmp(:,:) = sf(jp_cal)%fnow(:,:,1) * lfus               ! add the latent heat of iceberg melting
1014         qns_tot(:,:) = qns_tot(:,:) - ztmp(:,:)
1015         IF( lk_diaar5 )   CALL iom_put( 'hflx_cal_cea', ztmp + sf(jp_cal)%fnow(:,:,1) * zcptn(:,:) )   ! heat flux from calving
1016      ENDIF
1017
1018      !                                                      ! ========================= !
1019      SELECT CASE( TRIM( sn_rcv_qsr%cldes ) )                !      solar heat fluxes    !   (qsr)
1020      !                                                      ! ========================= !
1021      CASE( 'oce only' )
1022         qsr_tot(:,:  ) = MAX( 0._wp , sf(jp_qsroce)%fnow(:,:,1) )
1023      CASE( 'conservative' )
1024         qsr_tot(:,:  ) = sf(jp_qsrmix)%fnow(:,:,1)
1025         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
1026            qsr_ice(:,:,1:jpl) = sf(jp_qsrice)%fnow(:,:,1:jpl)
1027         ELSE
1028            ! Set all category values equal for the moment
1029            DO jl=1,jpl
1030               qsr_ice(:,:,jl) = sf(jp_qsrice)%fnow(:,:,1)
1031            ENDDO
1032         ENDIF
1033         qsr_tot(:,:  ) = sf(jp_qsrmix)%fnow(:,:,1)
1034         qsr_ice(:,:,1) = sf(jp_qsrice)%fnow(:,:,1)
1035      CASE( 'oce and ice' )
1036         qsr_tot(:,:  ) =  p_frld(:,:) * sf(jp_qsroce)%fnow(:,:,1)
1037         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
1038            DO jl=1,jpl
1039               qsr_tot(:,:   ) = qsr_tot(:,:) + a_i(:,:,jl) * sf(jp_qsrice)%fnow(:,:,jl)   
1040               qsr_ice(:,:,jl) = sf(jp_qsrice)%fnow(:,:,jl)
1041            ENDDO
1042         ELSE
1043            DO jl=1,jpl
1044               qsr_tot(:,:   ) = qsr_tot(:,:) + zicefr(:,:) * sf(jp_qsrice)%fnow(:,:,1)
1045               qsr_ice(:,:,jl) = sf(jp_qsrice)%fnow(:,:,1)
1046            ENDDO
1047         ENDIF
1048      CASE( 'mixed oce-ice' )
1049         qsr_tot(:,:  ) = sf(jp_qsrmix)%fnow(:,:,1)
1050! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
1051!       Create solar heat flux over ice using incoming solar heat flux and albedos
1052!       ( see OASIS3 user guide, 5th edition, p39 )
1053         qsr_ice(:,:,1) = sf(jp_qsrmix)%fnow(:,:,1) * ( 1.- palbi(:,:,1) )   &
1054            &            / (  1.- ( albedo_oce_mix(:,:  ) * p_frld(:,:)       &
1055            &                     + palbi         (:,:,1) * zicefr(:,:) ) )
1056      END SELECT
1057      IF( ln_dm2dc ) THEN   ! modify qsr to include the diurnal cycle
1058         qsr_tot(:,:  ) = sbc_dcy( qsr_tot(:,:  ) )
1059         DO jl=1,jpl
1060            qsr_ice(:,:,jl) = sbc_dcy( qsr_ice(:,:,jl) )
1061         ENDDO
1062      ENDIF
1063
1064      SELECT CASE( TRIM( sn_rcv_dqnsdt%cldes ) )
1065      CASE ('coupled')
1066         IF ( TRIM(sn_rcv_dqnsdt%clcat) == 'yes' ) THEN
1067            dqns_ice(:,:,1:jpl) = sf(jp_dqnsdt)%fnow(:,:,1:jpl)
1068         ELSE
1069            ! Set all category values equal for the moment
1070            DO jl=1,jpl
1071               dqns_ice(:,:,jl) = sf(jp_dqnsdt)%fnow(:,:,1)
1072            ENDDO
1073         ENDIF
1074      END SELECT
1075
1076      SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) )
1077      CASE ('coupled')
1078         topmelt(:,:,:)=sf(jp_topm)%fnow(:,:,:)
1079         botmelt(:,:,:)=sf(jp_botm)%fnow(:,:,:)
1080      END SELECT
1081
1082      CALL wrk_dealloc( jpi,jpj, zcptn, ztmp, zicefr )
1083      !
1084      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_ice_flx')
1085      !
1086   END SUBROUTINE sbc_cpl_ice_flx
1087   
1088   
1089   SUBROUTINE sbc_cpl_snd( kt )
1090      !!----------------------------------------------------------------------
1091      !!             ***  ROUTINE sbc_cpl_snd  ***
1092      !!
1093      !! ** Purpose :   provide the ocean-ice informations to the atmosphere
1094      !!
1095      !! ** Method  :   send to the atmosphere through a call to cpl_prism_snd
1096      !!              all the needed fields (as defined in sbc_cpl_init)
1097      !!----------------------------------------------------------------------
1098      INTEGER, INTENT(in) ::   kt
1099      !
1100      INTEGER ::   ji, jj, jl   ! dummy loop indices
1101      INTEGER ::   isec, info   ! local integer
1102      REAL(wp), POINTER, DIMENSION(:,:)   ::   zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1
1103      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmp3, ztmp4   
1104      !!----------------------------------------------------------------------
1105      !
1106      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_snd')
1107      !
1108      CALL wrk_alloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
1109      CALL wrk_alloc( jpi,jpj,jpl, ztmp3, ztmp4 )
1110
1111      isec = ( kt - nit000 ) * NINT(rdttra(1))        ! date of exchanges
1112
1113      zfr_l(:,:) = 1.- fr_i(:,:)
1114
1115      !                                                      ! ------------------------- !
1116      !                                                      !    Surface temperature    !   in Kelvin
1117      !                                                      ! ------------------------- !
1118      IF( ssnd(jps_toce)%laction .OR. ssnd(jps_tice)%laction .OR. ssnd(jps_tmix)%laction ) THEN
1119         SELECT CASE( sn_snd_temp%cldes)
1120         CASE( 'oce only'             )   ;   ztmp1(:,:) =   tsn(:,:,1,jp_tem) + rt0
1121         CASE( 'weighted oce and ice' )   ;   ztmp1(:,:) = ( tsn(:,:,1,jp_tem) + rt0 ) * zfr_l(:,:)   
1122            SELECT CASE( sn_snd_temp%clcat )
1123            CASE( 'yes' )   
1124               ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
1125            CASE( 'no' )
1126               ztmp3(:,:,:) = 0.0
1127               DO jl=1,jpl
1128                  ztmp3(:,:,1) = ztmp3(:,:,1) + tn_ice(:,:,jl) * a_i(:,:,jl)
1129               ENDDO
1130            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
1131            END SELECT
1132         CASE( 'mixed oce-ice'        )   
1133            ztmp1(:,:) = ( tsn(:,:,1,1) + rt0 ) * zfr_l(:,:) 
1134            DO jl=1,jpl
1135               ztmp1(:,:) = ztmp1(:,:) + tn_ice(:,:,jl) * a_i(:,:,jl)
1136            ENDDO
1137         CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%cldes' )
1138         END SELECT
1139         IF( ssnd(jps_toce)%laction )   CALL cpl_prism_snd( jps_toce, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
1140         IF( ssnd(jps_tice)%laction )   CALL cpl_prism_snd( jps_tice, isec, ztmp3, info )
1141         IF( ssnd(jps_tmix)%laction )   CALL cpl_prism_snd( jps_tmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
1142      ENDIF
1143      !
1144      !                                                      ! ------------------------- !
1145      !                                                      !           Albedo          !
1146      !                                                      ! ------------------------- !
1147      IF( ssnd(jps_albice)%laction ) THEN                         ! ice
1148         ztmp3(:,:,1:jpl) = alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
1149         CALL cpl_prism_snd( jps_albice, isec, ztmp3, info )
1150      ENDIF
1151      IF( ssnd(jps_albmix)%laction ) THEN                         ! mixed ice-ocean
1152         ztmp1(:,:) = albedo_oce_mix(:,:) * zfr_l(:,:)
1153         DO jl=1,jpl
1154            ztmp1(:,:) = ztmp1(:,:) + alb_ice(:,:,jl) * a_i(:,:,jl)
1155         ENDDO
1156         CALL cpl_prism_snd( jps_albmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
1157      ENDIF
1158      !                                                      ! ------------------------- !
1159      !                                                      !  Ice fraction & Thickness !
1160      !                                                      ! ------------------------- !
1161      ! Send ice fraction field
1162      IF( ssnd(jps_fice)%laction ) THEN
1163         SELECT CASE( sn_snd_thick%clcat )
1164         CASE( 'yes' )   ;   ztmp3(:,:,1:jpl) =  a_i(:,:,1:jpl)
1165         CASE( 'no'  )   ;   ztmp3(:,:,1    ) = fr_i(:,:      )
1166         CASE default    ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
1167         END SELECT
1168         CALL cpl_prism_snd( jps_fice, isec, ztmp3, info )
1169      ENDIF
1170
1171      ! Send ice and snow thickness field
1172      IF( ssnd(jps_hice)%laction .OR. ssnd(jps_hsnw)%laction ) THEN
1173         SELECT CASE( sn_snd_thick%cldes)
1174         CASE( 'none'                  )       ! nothing to do
1175         CASE( 'weighted ice and snow' )   
1176            SELECT CASE( sn_snd_thick%clcat )
1177            CASE( 'yes' )   
1178               ztmp3(:,:,1:jpl) =  ht_i(:,:,1:jpl) * a_i(:,:,1:jpl)
1179               ztmp4(:,:,1:jpl) =  ht_s(:,:,1:jpl) * a_i(:,:,1:jpl)
1180            CASE( 'no' )
1181               ztmp3(:,:,:) = 0.0   ;  ztmp4(:,:,:) = 0.0
1182               DO jl=1,jpl
1183                  ztmp3(:,:,1) = ztmp3(:,:,1) + ht_i(:,:,jl) * a_i(:,:,jl)
1184                  ztmp4(:,:,1) = ztmp4(:,:,1) + ht_s(:,:,jl) * a_i(:,:,jl)
1185               ENDDO
1186            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
1187            END SELECT
1188         CASE( 'ice and snow'         )   
1189            ztmp3(:,:,1:jpl) = ht_i(:,:,1:jpl)
1190            ztmp4(:,:,1:jpl) = ht_s(:,:,1:jpl)
1191         CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%cldes' )
1192         END SELECT
1193         IF( ssnd(jps_hice)%laction )   CALL cpl_prism_snd( jps_hice, isec, ztmp3, info )
1194         IF( ssnd(jps_hsnw)%laction )   CALL cpl_prism_snd( jps_hsnw, isec, ztmp4, info )
1195      ENDIF
1196      !
1197#if defined key_cpl_carbon_cycle
1198      !                                                      ! ------------------------- !
1199      !                                                      !  CO2 flux from PISCES     !
1200      !                                                      ! ------------------------- !
1201      IF( ssnd(jps_co2)%laction )   CALL cpl_prism_snd( jps_co2, isec, RESHAPE ( oce_co2, (/jpi,jpj,1/) ) , info )
1202      !
1203#endif
1204      !                                                      ! ------------------------- !
1205      IF( ssnd(jps_ocx1)%laction ) THEN                      !      Surface current      !
1206         !                                                   ! ------------------------- !
1207         !   
1208         !                                                  j+1   j     -----V---F
1209         ! surface velocity always sent from T point                     !       |
1210         !                                                        j      |   T   U
1211         !                                                               |       |
1212         !                                                   j    j-1   -I-------|
1213         !                                               (for I)         |       |
1214         !                                                              i-1  i   i
1215         !                                                               i      i+1 (for I)
1216         SELECT CASE( TRIM( sn_snd_crt%cldes ) )
1217         CASE( 'oce only'             )      ! C-grid ==> T
1218            DO jj = 2, jpjm1
1219               DO ji = fs_2, fs_jpim1   ! vector opt.
1220                  zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj  ,1) )
1221                  zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji  ,jj-1,1) ) 
1222               END DO
1223            END DO
1224         CASE( 'weighted oce and ice' )   
1225            SELECT CASE ( cp_ice_msh )
1226            CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
1227               DO jj = 2, jpjm1
1228                  DO ji = fs_2, fs_jpim1   ! vector opt.
1229                     zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj) 
1230                     zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)
1231                     zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
1232                     zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
1233                  END DO
1234               END DO
1235            CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
1236               DO jj = 2, jpjm1
1237                  DO ji = 2, jpim1   ! NO vector opt.
1238                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj) 
1239                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj) 
1240                     zitx1(ji,jj) = 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &
1241                        &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1242                     zity1(ji,jj) = 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &
1243                        &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1244                  END DO
1245               END DO
1246            CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
1247               DO jj = 2, jpjm1
1248                  DO ji = 2, jpim1   ! NO vector opt.
1249                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj) 
1250                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj) 
1251                     zitx1(ji,jj) = 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &
1252                        &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1253                     zity1(ji,jj) = 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &
1254                        &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1255                  END DO
1256               END DO
1257            END SELECT
1258            CALL lbc_lnk( zitx1, 'T', -1. )   ;   CALL lbc_lnk( zity1, 'T', -1. )
1259         CASE( 'mixed oce-ice'        )
1260            SELECT CASE ( cp_ice_msh )
1261            CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
1262               DO jj = 2, jpjm1
1263                  DO ji = fs_2, fs_jpim1   ! vector opt.
1264                     zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)   &
1265                        &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
1266                     zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)   &
1267                        &         + 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
1268                  END DO
1269               END DO
1270            CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
1271               DO jj = 2, jpjm1
1272                  DO ji = 2, jpim1   ! NO vector opt.
1273                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
1274                        &         + 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &
1275                        &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1276                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
1277                        &         + 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &
1278                        &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1279                  END DO
1280               END DO
1281            CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
1282               DO jj = 2, jpjm1
1283                  DO ji = 2, jpim1   ! NO vector opt.
1284                     zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
1285                        &         + 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &
1286                        &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
1287                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
1288                        &         + 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &
1289                        &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
1290                  END DO
1291               END DO
1292            END SELECT
1293         END SELECT
1294         CALL lbc_lnk( zotx1, ssnd(jps_ocx1)%clgrid, -1. )   ;   CALL lbc_lnk( zoty1, ssnd(jps_ocy1)%clgrid, -1. )
1295         !
1296         !
1297         IF( TRIM( sn_snd_crt%clvor ) == 'eastward-northward' ) THEN             ! Rotation of the components
1298            !                                                                     ! Ocean component
1299            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->e', ztmp1 )       ! 1st component
1300            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->n', ztmp2 )       ! 2nd component
1301            zotx1(:,:) = ztmp1(:,:)                                                   ! overwrite the components
1302            zoty1(:,:) = ztmp2(:,:)
1303            IF( ssnd(jps_ivx1)%laction ) THEN                                     ! Ice component
1304               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 )    ! 1st component
1305               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 )    ! 2nd component
1306               zitx1(:,:) = ztmp1(:,:)                                                ! overwrite the components
1307               zity1(:,:) = ztmp2(:,:)
1308            ENDIF
1309         ENDIF
1310         !
1311         ! spherical coordinates to cartesian -> 2 components to 3 components
1312         IF( TRIM( sn_snd_crt%clvref ) == 'cartesian' ) THEN
1313            ztmp1(:,:) = zotx1(:,:)                     ! ocean currents
1314            ztmp2(:,:) = zoty1(:,:)
1315            CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 )
1316            !
1317            IF( ssnd(jps_ivx1)%laction ) THEN           ! ice velocities
1318               ztmp1(:,:) = zitx1(:,:)
1319               ztmp1(:,:) = zity1(:,:)
1320               CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 )
1321            ENDIF
1322         ENDIF
1323         !
1324         IF( ssnd(jps_ocx1)%laction )   CALL cpl_prism_snd( jps_ocx1, isec, RESHAPE ( zotx1, (/jpi,jpj,1/) ), info )   ! ocean x current 1st grid
1325         IF( ssnd(jps_ocy1)%laction )   CALL cpl_prism_snd( jps_ocy1, isec, RESHAPE ( zoty1, (/jpi,jpj,1/) ), info )   ! ocean y current 1st grid
1326         IF( ssnd(jps_ocz1)%laction )   CALL cpl_prism_snd( jps_ocz1, isec, RESHAPE ( zotz1, (/jpi,jpj,1/) ), info )   ! ocean z current 1st grid
1327         !
1328         IF( ssnd(jps_ivx1)%laction )   CALL cpl_prism_snd( jps_ivx1, isec, RESHAPE ( zitx1, (/jpi,jpj,1/) ), info )   ! ice   x current 1st grid
1329         IF( ssnd(jps_ivy1)%laction )   CALL cpl_prism_snd( jps_ivy1, isec, RESHAPE ( zity1, (/jpi,jpj,1/) ), info )   ! ice   y current 1st grid
1330         IF( ssnd(jps_ivz1)%laction )   CALL cpl_prism_snd( jps_ivz1, isec, RESHAPE ( zitz1, (/jpi,jpj,1/) ), info )   ! ice   z current 1st grid
1331         !
1332      ENDIF
1333      !
1334      CALL wrk_dealloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
1335      CALL wrk_dealloc( jpi,jpj,jpl, ztmp3, ztmp4 )
1336      !
1337      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_snd')
1338      !
1339   END SUBROUTINE sbc_cpl_snd
1340   
1341#else
1342   !!----------------------------------------------------------------------
1343   !!   Dummy module                                            NO coupling
1344   !!----------------------------------------------------------------------
1345   USE par_kind        ! kind definition
1346CONTAINS
1347   SUBROUTINE sbc_cpl_snd( kt )
1348      WRITE(*,*) 'sbc_cpl_snd: You should not have seen this print! error?', kt
1349   END SUBROUTINE sbc_cpl_snd
1350   !
1351   SUBROUTINE sbc_cpl_rcv( kt, k_fsbc, k_ice )     
1352      WRITE(*,*) 'sbc_cpl_snd: You should not have seen this print! error?', kt, k_fsbc, k_ice
1353   END SUBROUTINE sbc_cpl_rcv
1354   !
1355   SUBROUTINE sbc_cpl_ice_tau( p_taui, p_tauj )     
1356      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_taui   ! i- & j-components of atmos-ice stress [N/m2]
1357      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_tauj   ! at I-point (B-grid) or U & V-point (C-grid)
1358      p_taui(:,:) = 0.   ;   p_tauj(:,:) = 0. ! stupid definition to avoid warning message when compiling...
1359      WRITE(*,*) 'sbc_cpl_snd: You should not have seen this print! error?'
1360   END SUBROUTINE sbc_cpl_ice_tau
1361   !
1362   SUBROUTINE sbc_cpl_ice_flx( p_frld , palbi   , psst    , pist  )
1363      REAL(wp), INTENT(in   ), DIMENSION(:,:  ) ::   p_frld     ! lead fraction                [0 to 1]
1364      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   palbi   ! ice albedo
1365      REAL(wp), INTENT(in   ), DIMENSION(:,:  ), OPTIONAL ::   psst    ! sea surface temperature      [Celcius]
1366      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   pist    ! ice surface temperature      [Kelvin]
1367      WRITE(*,*) 'sbc_cpl_snd: You should not have seen this print! error?', p_frld(1,1), palbi(1,1,1), psst(1,1), pist(1,1,1) 
1368   END SUBROUTINE sbc_cpl_ice_flx
1369   
1370#endif
1371
1372   !!======================================================================
1373END MODULE sbccpl
Note: See TracBrowser for help on using the repository browser.