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.
ocesbc.F90 in trunk/NEMO/OPA_SRC/SBC – NEMO

source: trunk/NEMO/OPA_SRC/SBC/ocesbc.F90 @ 700

Last change on this file since 700 was 700, checked in by smasson, 17 years ago

sbc(1): bugfix on wind stress computation #2

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 36.7 KB
Line 
1MODULE ocesbc
2   !!======================================================================
3   !!                     ***  MODULE  ocesbc  ***
4   !!                     Ocean surface boundary conditions
5   !!======================================================================
6
7   !!----------------------------------------------------------------------
8   !!   oce_sbc     : ???
9   !!   oce_sbc_dmp : ???
10   !!----------------------------------------------------------------------
11   !! * Modules used
12   USE oce            ! dynamics and tracers variables
13   USE dom_oce        ! ocean space domain variables
14   USE ice_oce        ! sea-ice variable
15   USE blk_oce        ! bulk variables
16   USE flx_oce        ! sea-ice/ocean forcings variables
17   USE phycst         ! Define parameters for the routines
18   USE taumod         ! surface stress forcing
19   USE flxmod         ! thermohaline fluxes
20   USE flxrnf         ! runoffs forcing
21   USE tradmp         ! damping salinity trend
22   USE dtatem         ! ocean temperature data
23   USE dtasal         ! ocean salinity data
24   USE ocfzpt         ! surface ocean freezing point
25   USE lbclnk         ! ocean lateral boundary condition
26   USE lib_mpp        ! distribued memory computing library
27   USE in_out_manager ! I/O manager
28   USE prtctl         ! Print control
29   USE diurnalcycle   !
30
31   IMPLICIT NONE
32   PRIVATE
33
34   !! * Accessibility
35   PUBLIC oce_sbc    ! routine called by step
36
37   !! * Shared module variables
38   REAL(wp), PUBLIC ::   &  !:
39      aplus, aminus,     &  !:
40      empold = 0.e0,     &  !: current year freshwater budget correction
41      area_g                !: surface og the global ocean
42   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !:
43!CT      qt  ,         &  !: total surface heat flux (w/m2)
44!CT      qsr ,         &  !: solar radiation (w/m2)
45      qla ,         &  !: Latent heat flux (W/m2)
46      qlw ,         &  !: Long Wave Heat Flux (W/m2)
47      qsb ,         &  !: Sensible Heat Flux (W/m2)
48      emp ,         &  !: evaporation minus precipitation (kg/m2/s = mm/s)
49      emps,         &  !: evaporation - precipitation (free surface)
50      qrp ,         &  !: heat flux damping (w/m2)
51      erp              !: evaporation damping (kg/m2/s = mm/s)
52#if ! defined key_dynspg_rl
53   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !:
54      dmp              !: internal dampind term
55   REAL(wp), DIMENSION(jpi,jpj) ::   &  !:
56      e1e2_i    ! area of the interior domain (e1t*e2t*tmask_i)
57#endif
58
59#  include "domzgr_substitute.h90"
60#  include "vectopt_loop_substitute.h90"
61
62   !!----------------------------------------------------------------------
63   !!  OPA 9.0 , LOCEAN-IPSL (2005)
64   !! $Header$
65   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
66   !!----------------------------------------------------------------------
67CONTAINS
68
69#if defined key_ice_lim
70   !!----------------------------------------------------------------------
71   !!   'key_ice_lim' :                                   LIM sea-ice model
72   !!----------------------------------------------------------------------
73# if defined key_coupled
74      !!----------------------------------------------------------------------
75      !!   'key_coupled' :                            Coupled Ocean/Atmosphere
76      !!----------------------------------------------------------------------
77
78   SUBROUTINE oce_sbc( kt )
79      !!---------------------------------------------------------------------
80      !!                 ***  ROUTINE oce_sbc  ***
81      !!                   
82      !! ** Purpose :   Ocean surface boundaries conditions with
83      !!        Louvain la Neuve Sea Ice Model in coupled mode
84      !!
85      !! History :
86      !!   1.0  !  00-10  (O. Marti)  Original code
87      !!   2.0  !  02-12  (G. Madec)  F90: Free form and module
88      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
89      !!----------------------------------------------------------------------
90      !! * Arguments
91      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
92
93      !! * Local declarations
94      INTEGER ::   ji, jj                   ! dummy loop indices
95      REAL(wp) ::   ztx, ztaux_ice, zty, ztauy_ice
96      REAL(wp) ::   ztdta, ztgel, zqrp, zm_emp
97      !!----------------------------------------------------------------------
98
99      IF( kt == nit000 ) THEN
100         e1e2_i(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:)
101         area_g = SUM( e1e2_i(:,:) )
102         IF( lk_mpp )   CALL  mpp_sum( area_g )   ! sum over the global domain
103      ENDIF
104
105      ! 1. initialization to zero at kt = nit000
106      ! ---------------------------------------
107
108      IF( MOD( kt-1, nfice ) == 0 ) THEN 
109
110         ! non solar heat flux + solar flux + restoring
111         qt (:,:) = fnsolar(:,:) + fsolar(:,:)
112
113         ! solar flux
114         qsr(:,:) = fsolar(:,:)
115
116#if ! defined key_dynspg_rl 
117         ! total concentration/dilution effect (use on SSS)
118         emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + calving(:,:)
119         ! total volume flux (use on sea-surface height)
120         emp (:,:) = fmass(:,:)              + runoff(:,:) + calving(:,:)
121
122         ! Compute the mean emp(:,:)
123         zm_emp =  SUM( emp(:,:)*e1e2_i(:,:) ) / area_g
124         IF( lk_mpp )   CALL  mpp_sum( zm_emp )   ! sum over the global domain
125
126         ! Correct emp() & emps() in substracting the emp() mean
127         emps(:,:) = emps(:,:) - zm_emp
128         emp (:,:) = emp (:,:) - zm_emp
129#else
130         ! Rigid-lid (emp=emps=E-P-R+Erp)
131         ! freshwater flux
132         emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + calving(:,:)
133         emp (:,:) = emps(:,:)
134#endif
135
136         DO jj = 1, jpjm1
137            DO ji = 1, fs_jpim1   ! vertor opt.
138               ztx         = 0.5 * ( freeze(ji  ,jj) + freeze(ji+1,jj  ) )   !   ice cover        at U-point
139               ztaux_ice   = 0.5 * ( ftaux (ji+1,jj) + ftaux (ji+1,jj+1) )   !   ice/ocean stress at U-point
140               taux(ji,jj) = (1.0-ztx) * taux_oce(ji,jj) + ztx * ztaux_ice
141
142               zty         = 0.5 * ( freeze(ji,jj  ) + freeze(ji  ,jj+1) )   !   ice cover        at V-point
143               ztauy_ice   = 0.5 * ( ftauy (ji,jj+1) + ftauy (ji+1,jj+1) )   !   ice/ocean stress at V-point
144               tauy(ji,jj) = (1.0-zty) * tauy_oce(ji,jj) + zty * ztauy_ice
145            END DO
146         END DO
147         CALL lbc_lnk( taux, 'U', -1. )
148         CALL lbc_lnk( tauy, 'V', -1. )   
149
150         ! Re-initialization of fluxes
151         sst_io(:,:) = 0.e0
152         sss_io(:,:) = 0.e0
153         u_io  (:,:) = 0.e0
154         v_io  (:,:) = 0.e0
155         gtaux (:,:) = 0.e0
156         gtauy (:,:) = 0.e0
157
158      ENDIF
159
160   END SUBROUTINE oce_sbc
161
162# elif defined key_flx_bulk_monthly || defined key_flx_bulk_daily || defined key_flx_core
163      !!----------------------------------------------------------------------
164      !!   'key_ice_lim'                              with  LIM sea-ice model
165      !!----------------------------------------------------------------------
166
167   SUBROUTINE oce_sbc( kt )
168      !!---------------------------------------------------------------------
169      !!                   ***  ROUTINE oce_sbc  ***
170      !!                   
171      !! ** Purpose : - Ocean surface boundary conditions with LIM sea-ice
172      !!        model in forced mode using bulk formulea
173      !!
174      !! History :
175      !!   1.0  !  99-11  (M. Imbard)  Original code
176      !!        !  01-03  (D. Ludicone, E. Durand, G. Madec) free surf.
177      !!   2.0  !  02-09  (G. Madec, C. Ethe)  F90: Free form and module
178      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
179      !!----------------------------------------------------------------------
180      !! * arguments
181      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
182
183      !! * Local declarations
184      INTEGER  ::   ji, jj                   ! dummy loop indices
185      REAL(wp) ::   ztx, ztaux, zty, ztauy, zm_emp
186      !!----------------------------------------------------------------------
187
188      ! 1. initialization to zero at kt = nit000
189      ! ---------------------------------------
190
191      IF( kt == nit000 ) THEN     
192         qsr    (:,:) = 0.e0
193         qt     (:,:) = 0.e0
194         qla    (:,:) = 0.e0
195         qlw    (:,:) = 0.e0
196         qsb    (:,:) = 0.e0
197         qrp    (:,:) = 0.e0
198         emp    (:,:) = 0.e0
199         emps   (:,:) = 0.e0
200         erp    (:,:) = 0.e0
201#if ! defined key_dynspg_rl
202         dmp    (:,:) = 0.e0
203#endif
204         e1e2_i(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:)
205         area_g = SUM( e1e2_i(:,:) )
206
207         IF( lk_mpp )   CALL  mpp_sum( area_g )   ! sum over the global domain
208      ENDIF
209#if defined key_flx_core
210      qla(:,:) = (1 - freeze(:,:))*qla_oce + freeze(:,:)*qla_ice(:,:)
211      qsb(:,:) = (1 - freeze(:,:))*qsb_oce + freeze(:,:)*qsb_ice(:,:)
212      qlw(:,:) = (1 - freeze(:,:))*qlw_oce + freeze(:,:)*qlw_ice(:,:)
213#endif
214      IF( MOD( kt-1, nfice ) == 0 ) THEN
215
216         CALL oce_sbc_dmp       ! Computation of internal and evaporation damping terms       
217
218         ! Surface Ocean fluxes
219         ! ====================
220
221         ! Surface heat flux (W/m2)
222         ! -----------------
223
224         qt  (:,:) = fnsolar(:,:) + fsolar(:,:)     ! non solar heat flux + solar flux
225         qsr (:,:) = fsolar(:,:)                     ! solar flux
226
227#if ! defined key_dynspg_rl     
228         ! total concentration/dilution effect (use on SSS)
229         emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + erp(:,:) + empold
230
231         ! total volume flux (use on sea-surface height)
232         emp (:,:) = fmass(:,:) -   dmp(:,:) + runoff(:,:) + erp(:,:) + empold     
233
234         ! Compute the mean emp(:,:)
235         zm_emp =  SUM( emp(:,:)*e1e2_i(:,:) ) / area_g
236
237         IF( lk_mpp )   CALL  mpp_sum( zm_emp )   ! sum over the global domain
238
239         ! Correct emp() & emps() in substracting the emp() mean
240         emps(:,:) = emps(:,:) - zm_emp
241         emp (:,:) = emp (:,:) - zm_emp
242#else
243         ! Rigid-lid (emp=emps=E-P-R+Erp)
244         emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + erp(:,:)     ! freshwater flux
245         emp (:,:) = emps(:,:)
246
247#endif
248
249         ! Surface stress
250         ! --------------
251
252         ! update the stress beloww sea-ice area
253         DO jj = 1, jpjm1
254            DO ji = 1, fs_jpim1   ! vertor opt.
255               ztx         = MAX( freezn(ji,jj), freezn(ji,jj+1) )   ! ice/ocean indicator at U- and V-points
256               zty         = MAX( freezn(ji,jj), freezn(ji+1,jj) )
257               ztaux       = 0.5 *( ftaux(ji+1,jj) + ftaux(ji+1,jj+1) ) ! ice-ocean stress at U- and V-points
258               ztauy       = 0.5 *( ftauy(ji,jj+1) + ftauy(ji+1,jj+1) )
259               taux(ji,jj) = (1.-ztx) * taux(ji,jj) + ztx * ztaux    ! stress at the ocean surface
260               tauy(ji,jj) = (1.-zty) * tauy(ji,jj) + zty * ztauy
261            END DO
262         END DO
263
264         ! boundary condition on the stress (taux,tauy)
265         CALL lbc_lnk( taux, 'U', -1. )
266         CALL lbc_lnk( tauy, 'V', -1. )
267
268         ! Re-initialization of fluxes
269         sst_io(:,:) = 0.e0
270         sss_io(:,:) = 0.e0
271         u_io  (:,:) = 0.e0
272         v_io  (:,:) = 0.e0
273
274      ENDIF
275
276      IF( ln_dcy )   CALL diurnal_cycle( kt )
277
278   END SUBROUTINE oce_sbc
279
280# else
281      !!----------------------------------------------------------------------
282      !!   Error option               LIM sea-ice model requires bulk formulea
283      !!----------------------------------------------------------------------
284      This line forced a compilation error
285# endif
286
287#else
288   !!----------------------------------------------------------------------
289   !!   Default option                                 NO LIM sea-ice model
290   !!----------------------------------------------------------------------
291# if defined key_coupled
292      !!----------------------------------------------------------------------
293      !!   'key_coupled' :                            Coupled Ocean/Atmosphere
294      !!----------------------------------------------------------------------
295
296   SUBROUTINE oce_sbc( kt )
297      !!---------------------------------------------------------------------
298      !!                 ***  ROUTINE oce_sbc  ***
299      !!                   
300      !! ** Purpose :   Ocean surface boundaries conditions in
301      !!                coupled ocean/atmosphere case without sea-ice
302      !!
303      !! History :
304      !!   1.0  !  00-10  (O. Marti)  Original code
305      !!   2.0  !  02-12  (G. Madec)  F90: Free form and module
306      !!----------------------------------------------------------------------
307      !! * Arguments
308      INTEGER, INTENT( in  ) ::   kt   ! ocean time step index
309
310      !! * Local declarations
311      INTEGER  ::   ji, jj, jf         ! dummy loop indices
312      REAL(wp) ::   ztgel,          &  ! temporary scalars
313         zice, zhemis, zqrp, zqri,  &  !    "         "
314         zq, zqi, zerp, ze, zei, zro   !    "         "
315      !!----------------------------------------------------------------------
316
317      ! Compute fluxes
318      ! --------------
319
320      DO jj = 1, jpj
321         DO ji = 1, jpi
322
323            ztgel = fzptn(ji,jj)   ! local freezing temperature
324
325            ! opa model ice freeze()
326
327            zice = tmask(ji,jj,1)
328            IF( tn(ji,jj,1) >=  ztgel ) zice = 0.
329            freeze(ji,jj) = zice
330
331            ! hemisphere indicator (=1 north, =-1 south)
332
333            zhemis = float(isign(1, mjg(jj)-(jpjglo/2+1)))
334
335            ! a) net downward radiative flux qsr()
336            ! - AGCM qsrc if no ice
337            ! - zero under ice (zice=1)
338
339            qsr(ji,jj) = (1.-zice)*qsrc(ji,jj)*tmask(ji,jj,1)
340
341            ! b) heat flux damping term qrp()
342            ! - no damping if no  ice      (zice=0)
343            ! - gamma*min(0,t-tgel) if ice (zice=1)
344
345            zqrp = 0.
346            zqri = dqdt0*MIN( 0., tb(ji,jj,1)-ztgel )
347            qrp(ji,jj) = ( ( 1. - zice ) * zqrp + zice * zqri ) * tmask(ji,jj,1)
348
349
350            ! c) net downward heat flux q() = q0 + qrp()
351            ! for q0
352            ! - AGCM qc if no  ice (zice=0)
353            ! - -2 watt/m2 (arctic) or -4 watt/m2 (antarctic) if ice (zice=1)
354            zq  = qc(ji,jj)
355            zqi = -3. + zhemis
356            qt(ji,jj) = ( (1.-zice) * zq + zice * zqi ) * tmask(ji,jj,1) + qrp(ji,jj)
357
358            ! d) water flux damping term erp()
359            ! - no damping
360            zerp = 0.
361            erp(ji,jj) = zerp
362
363            ! e) net upward water flux e() = eo + runoff() + erp()
364            ! for e0
365            ! - AGCM if no ice (zice=0)
366            ! - 1.mm/day if climatological and opa ice (zice=1)
367            ze  = ec(ji,jj)
368            zei = 1./rday
369            zro = runoff(ji,jj)
370            emp(ji,jj) = ( ( 1. - zice ) *  ze + zice * zei + zro ) * tmask(ji,jj,1) + erp(ji,jj)
371
372            ! f) net upward water flux for the salinity surface
373            !    boundary condition
374            emps(:,:) = emp(:,:)
375
376         END DO
377      END DO
378
379   END SUBROUTINE oce_sbc
380# elif defined key_flx_bulk_monthly || defined key_flx_bulk_daily || defined key_flx_forced_daily || defined key_flx_core
381      !!-------------------------------------------------------------------------
382      !!   'key_flx_bulk_monthly' or 'key_flx_bulk_daily' or  core   bulk formulea
383      !!   'key_flx_forced_daily'                                or no bulk case
384      !!-------------------------------------------------------------------------
385
386   SUBROUTINE oce_sbc( kt )
387      !!---------------------------------------------------------------------
388      !!                   ***  ROUTINE oce_sbc  ***
389      !!                   
390      !! ** Purpose :   Ocean surface boundary conditions in forced mode
391      !!      using either flux or bulk formulation.
392      !!
393      !! History :
394      !!   1.0  !  99-11  (M. Imbard)  Original code
395      !!        !  01-03  (D. Ludicone, E. Durand, G. Madec) free surf.
396      !!   2.0  !  02-09  (G. Madec, C. Ethe)  F90: Free form and module
397      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
398      !!----------------------------------------------------------------------
399      !! * Modules used
400      USE daymod                 ! calendar
401#if ! defined key_dtasst
402      USE dtasst, ONLY : rclice  ! sea surface temperature data
403#endif
404#if defined key_flx_bulk_monthly || defined key_flx_bulk_daily || defined key_flx_core
405      USE blk_oce                ! bulk variables
406#endif
407#if defined key_flx_forced_daily
408      USE flx_oce                ! sea-ice/ocean forcings variables
409#endif
410
411      !! * arguments
412      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
413
414      !! * local declarations
415      INTEGER ::   ji, jj        ! dummy loop arguments
416      INTEGER ::   i15, ifreq             !     
417      REAL(wp) ::  zxy
418      REAL(wp) ::  zsice, zqri, zqrp, ztdta, zqrj
419      REAL(wp) ::  zq, zqi, zhemis
420      REAL(wp), DIMENSION(jpi,jpj) :: zeri, zerps, ziclim
421      REAL(wp), DIMENSION(jpi,jpj) :: zqt, zqsr, zemp 
422      !!----------------------------------------------------------------------
423 
424
425#if defined key_flx_core
426      CALL ctl_stop( 'flxcore and no ice model not tested yet' )
427#endif
428
429      ! 1. initialization to zero at kt = nit000
430      ! ---------------------------------------
431
432      IF( kt == nit000 ) THEN     
433         qsr    (:,:) = 0.e0
434         freeze (:,:) = 0.e0
435         qt     (:,:) = 0.e0
436         qla    (:,:) = 0.e0
437         qlw    (:,:) = 0.e0
438         qsb    (:,:) = 0.e0
439         qrp    (:,:) = 0.e0
440         emp    (:,:) = 0.e0
441         emps   (:,:) = 0.e0
442         erp    (:,:) = 0.e0
443#if ! defined key_dynspg_rl
444         dmp    (:,:) = 0.e0
445#endif
446      ENDIF
447
448#if defined key_flx_bulk_monthly || defined key_flx_bulk_daily || defined key_flx_core
449      ifreq      = nfbulk
450      zqt (:,:)  = qsr_oce(:,:) + qnsr_oce(:,:)
451      zqsr(:,:)  = qsr_oce(:,:)
452      zemp(:,:)  = evap(:,:) - tprecip(:,:)
453#endif 
454
455#if defined key_flx_forced_daily
456      ifreq      = 1
457      zqt (:,:)  = p_qt (:,:)
458      zqsr(:,:)  = p_qsr(:,:)
459      zemp(:,:)  = p_emp(:,:)
460#endif 
461#if defined key_flx_core
462      qla(:,:) = (1 - freeze(:,:))*qla_oce + freeze(:,:)*qla_ice(:,:)
463      qsb(:,:) = (1 - freeze(:,:))*qsb_oce + freeze(:,:)*qsb_ice(:,:)
464      qlw(:,:) = (1 - freeze(:,:))*qlw_oce + freeze(:,:)*qlw_ice(:,:)
465#endif
466      IF( MOD( kt-1, ifreq) == 0 ) THEN
467         ! Computation of internal and evaporation damping terms       
468         CALL oce_sbc_dmp
469
470         zsice = - 0.04 / 0.8    ! ratio of isohaline compressibility over isotherme compressibility
471                                 ! ( d rho / dt ) / ( d rho / ds )      ( s = 34, t = -1.8 )
472         ! Flux computation
473         DO jj = 1, jpj
474            DO ji = 1, jpi     
475               ! climatological ice
476               ziclim(ji,jj) = FLOAT( NINT( rclice(ji,jj,1) ) )
477
478               ! avoid surfreezing point           
479               tn(ji,jj,1) = MAX( tn(ji,jj,1), fzptn(ji,jj) )
480
481               ! hemisphere indicator (=1 north, =-1 south)           
482               zhemis = FLOAT( isign(1, mjg(jj) - (jpjdta/2+1) ) )
483
484               ! restoring temperature (ztdta >= to local freezing temperature)           
485#if defined key_dtasst
486               ztdta = MAX( sst(ji,jj),    fzptn(ji,jj) )
487#else
488               ztdta = MAX( t_dta(ji,jj,1), fzptn(ji,jj) )
489#endif
490
491               ! a) net downward radiative flux qsr()           
492               qsr(ji,jj) = (1.-ziclim(ji,jj)) * zqsr(ji,jj) * tmask(ji,jj,1)
493
494               ! b) heat flux damping term qrp()
495               ! - gamma*(t-tlevitus) if no  climatological ice (ziclim=0)
496               ! - gamma*(t-(tgel-1.))  if climatological ice and no opa ice   (ziclim=1 zicopa=0)
497               ! - gamma*min(0,t-tgel) if climatological and opa ice (ziclim=1 zicopa=1)
498
499               zqri = dqdt0 * ( tb(ji,jj,1) - ( fzptn(ji,jj) - 1.) )
500               zqrj = dqdt0 * MIN( 0., tb(ji,jj,1) - fzptn(ji,jj) )
501
502               qrp(ji,jj) =  ( ziclim(ji,jj) * ( (1 - freeze(ji,jj)) * zqri    &
503                 &                                  + freeze(ji,jj)  * zqrj ) ) * tmask(ji,jj,1)
504
505#if ! defined key_flx_bulk_monthly && ! defined key_flx_bulk_daily && ! defined key_flx_core
506               zqrp = dqdt0 * ( tb(ji,jj,1) - ztdta )
507               qrp(ji,jj) = qrp(ji,jj) + (1. - ziclim(ji,jj)) * zqrp
508# endif
509
510               ! c) net downward heat flux q() = q0 + qrp()
511               ! for q0
512               ! - ECMWF fluxes if no climatological ice      (ziclim=0)
513               ! - qrp if climatological ice and no opa ice   (ziclim=1 zicopa=0)
514               ! - -2 watt/m2 (arctic) or -4 watt/m2 (antarctic) if climatological and opa ice
515               !                                              (ziclim=1 zicopa=1)
516               zq  = zqt(ji,jj)
517               zqi = -3. + zhemis
518               qt (ji,jj) = ( (1.-ziclim(ji,jj)) * zq   &
519                  +ziclim(ji,jj)  * freeze(ji,jj) * zqi )   &
520                  * tmask(ji,jj,1)   &
521                  + qrp(ji,jj)
522
523            END DO
524         END DO
525
526#if ! defined key_dynspg_rl
527         ! Free-surface
528
529         ! Water flux for zero buoyancy flux if no opa ice and ice clim
530         zeri(:,:) = -zsice * qrp(:,:) * ro0cpr * rauw / 34.0
531         zerps(:,:) = ziclim(:,:) * ( (1-freeze(:,:)) * zeri(:,:) )
532
533         ! Contribution to sea level:
534         ! net upward water flux emp() = e-p + runoff() + erp() + dmp + empold
535         emp (:,:) = zemp(:,:)     &   ! e-p data
536            &      + runoff(:,:)   &   ! runoff data
537            &      + erp(:,:)      &   ! restoring term to SSS data
538            &      + dmp(:,:)      &   ! freshwater flux associated with internal damping
539            &      + empold            ! domain averaged annual mean correction
540
541         ! Contribution to salinity:
542         ! net upward water flux emps() = e-p + runoff() + erp() + zerps + empold
543         emps(:,:) = zemp(:,:)     &
544            &      + runoff(:,:)   &
545            &      + erp(:,:)      &
546            &      + zerps(:,:)    &
547            &      + empold 
548#else
549         ! Rigid-lid (emp=emps=E-P-R+Erp)
550         ! freshwater flux
551         zeri(:,:)  = -zsice * qrp(:,:) * ro0cpr * rauw / 34.0
552         zerps(:,:) = ziclim(:,:) * ( (1-freeze(:,:)) * zeri(:,:) )
553         emps (:,:) = zemp(:,:)     &
554            &       + runoff(:,:)   &
555            &       + erp(:,:)      &
556            &       + zerps(:,:)
557         emp (:,:) = emps(:,:)
558#endif 
559
560         ! Boundary condition on emp for free surface option
561         ! -------------------------------------------------
562         CALL lbc_lnk( emp, 'T', 1. )
563
564      ENDIF
565
566   END SUBROUTINE oce_sbc
567
568# else
569      !!----------------------------------------------------------------------
570      !!   Default option :                                 Analytical forcing
571      !!----------------------------------------------------------------------
572
573   SUBROUTINE oce_sbc( kt )
574      !!---------------------------------------------------------------------
575      !!                    ***  ROUTINE oce_sbc  ***
576      !!             
577      !! ** Purpose :   provide the thermohaline fluxes (heat and freshwater)
578      !!                to the ocean at each time step.
579      !!
580      !! ** Method  :   Constant surface fluxes (read in namelist (namflx))
581      !!
582      !! ** Action  : - qt, qsr, emp, emps, qrp, erp
583      !!
584      !! History :
585      !!        !  91-03  ()  Original code
586      !!   8.5  !  02-09  (G. Madec)  F90: Free form and module
587      !!   9.0  !  04-05  (A. Koch-Larrouy) Add Gyre configuration
588      !!----------------------------------------------------------------------
589      !! * Modules used
590      USE flxrnf                       ! ocean runoffs
591      USE daymod, ONLY : nyear         ! calendar
592      USE dtasss                       ! sea surface salinity data
593
594      !! * arguments
595      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
596
597      !! * local declarations
598      REAL(wp) ::                   & !!! surface fluxes namelist (namflx)
599         q0   = 0.e0,               &  ! net heat flux
600         qsr0 = 0.e0,               &  ! solar heat flux
601         emp0 = 0.e0                   ! net freshwater flux
602      REAL(wp) ::   zsrp,           &
603         zemp_S, zemp_N, zemp_sais, &
604         zTstar, zcos_sais1, zconv, &
605         zcos_sais2
606      REAL(wp) ::           &
607         zsumemp,           &          ! tampon used for the emp sum
608         zsurf,             &          ! tampon used for the domain sum
609         ztime,             &          ! time in hour
610         ztimemax1, ztimemin1, &       ! 21th june,   and 21th december if date0 = 1st january
611         ztimemax2, ztimemin2          ! 21th august, and 21th february if date0 = 1st january
612      REAL(wp), DIMENSION(jpi,jpj) :: t_star
613      INTEGER  ::   ji, jj             ! dummy loop indices
614
615      INTEGER  ::           &
616         zyear0,            &          ! initial year
617         zmonth0,           &          ! initial month
618         zday0,             &          ! initial day
619         zday_year0                    ! initial day since january 1st
620
621      NAMELIST/namflx/ q0, qsr0, emp0
622      !!---------------------------------------------------------------------
623
624      !same temperature, E-P as in HAZELEGER 2000
625
626      IF( cp_cfg == 'gyre' ) THEN
627
628          zyear0  = ndate0 / 10000                              ! initial year
629          zmonth0 = ( ndate0 - zyear0 * 10000 ) / 100           ! initial month
630          zday0   = ndate0 - zyear0 * 10000 - zmonth0 * 100     ! initial day betwen 1 and 30
631          zday_year0= (zmonth0-1)*30.+zday0                     ! initial day betwen 1 and 360
632
633          ! current day (in hours) since january the 1st of the current year
634          ztime = FLOAT( kt ) * rdt / (rmmss * rhhmm)  &    !  total incrementation (in hours)
635             &     - (nyear  - 1) * rjjhh * raajj           !  minus years since beginning of experiment (in hours)
636
637          ! 21th june at 24h  in hours
638          ztimemax1 = ((5.*30.)+21.)* 24.
639          ! 21th december in hours
640          ! rjjhh * raajj / 4 = 1 seasonal cycle in hours
641          ztimemin1 = ztimemax1 + rjjhh * raajj / 2         
642          ! 21th july at 24h in hours
643          ztimemax2 = ((6.*30.)+21.)* 24.
644          ! 21th january day in hours
645          ! rjjhh * raajj / 4 = 1 seasonal cycle in hours
646          ztimemin2 = ztimemax2 - rjjhh * raajj / 2         
647
648          ! amplitudes
649          zemp_S   = 0.7       ! intensity of COS in the South
650          zemp_N   = 0.8       ! intensity of COS in the North
651          zemp_sais= 0.1
652          zTstar   = 28.3      ! intemsity from 28.3 a -5 deg
653
654          ! 1/2 period between 21th June and 21th December
655          zcos_sais1 = COS( (ztime - ztimemax1) / (ztimemin1 - ztimemax1) * rpi )   
656
657          ! 1/2 period between 21th July and 21th January
658          zcos_sais2 = COS( (ztime - ztimemax2) / (ztimemax2 - ztimemin2) * rpi )
659
660          zconv = 3.16e-5      ! convert 1m/yr->3.16e-5mm/s
661          DO jj = 1, jpj
662             DO ji = 1, jpi
663                ! domain from 15 deg to 50 deg between 27 and 28  degC at 15N, -3
664                ! and 13 degC at 50N 53.5 + or - 11 = 1/4 period :
665                ! 64.5 in summer, 42.5 in winter
666                t_star (ji,jj) = zTstar * ( 1 + 1. / 50. * zcos_sais2 )   &                 
667                   &             * COS( rpi * (gphit(ji,jj) - 5.)        &                 
668                   &                        / (53.5 * ( 1 + 11 / 53.5 * zcos_sais2 ) * 2.) ) 
669                qt  (ji,jj) = dqdt0 * ( tb(ji,jj,1) - t_star(ji,jj) )
670                IF( gphit(ji,jj) >= 14.845 .AND. 37.2 >= gphit(ji,jj)) THEN
671                   ! zero at 37.8 deg, max at 24.6 deg
672                   emp  (ji,jj) =   zemp_S * zconv   &
673                      &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (24.6 - 37.2) )  & 
674                      &         * ( 1 - zemp_sais / zemp_S * zcos_sais1)
675                   emps (ji,jj) = emp  (ji,jj)
676                ELSE
677                   ! zero at 37.8 deg, max at 46.8 deg
678                   emp (ji,jj) =  - zemp_N * zconv   &
679                      &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (46.8 - 37.2) )  & 
680                      &         * ( 1 - zemp_sais / zemp_N * zcos_sais1 )
681                   emps (ji,jj) = emp  (ji,jj)
682                ENDIF
683                ! 23.5 deg : tropics
684                qsr (ji,jj) =  230 * COS( 3.1415 * ( gphit(ji,jj) - 23.5 * zcos_sais1 ) / ( 0.9 * 180 ) )   
685             END DO
686          END DO
687
688          ! Compute the emp flux such as its integration on the whole domain at each time is zero
689          IF( nbench /= 1 .AND. nbit_cmp /= 1 ) THEN
690             zsumemp = 0.e0   ;   zsurf = 0.e0
691             DO jj = 1, jpj
692                DO ji = 1, jpi
693                   zsumemp = zsumemp + emp(ji,jj) * tmask(ji,jj,1) * tmask_i(ji,jj)
694                   zsurf =  zsurf + tmask(ji,jj,1) * tmask_i(ji,jj)
695                END DO
696             END DO
697
698             IF( lk_mpp )   CALL mpp_sum( zsumemp  )       ! sum over the global domain
699             IF( lk_mpp )   CALL mpp_sum( zsurf    )       ! sum over the global domain
700
701             ! Default GYRE configuration
702             zsumemp = zsumemp / zsurf
703          ELSE
704             ! Benchmark GYRE configuration (to allow the bit to bit comparison between Mpp/Mono case)
705             zsumemp = 0.e0   ;    zsurf = 0.e0
706          ENDIF
707
708          IF( lk_dtasss ) THEN           ! Sea surface salinity damping
709             ! deds0 from (mm/day) => (mm/s)
710             zsrp = deds0 / rday 
711             erp(:,:) = zsrp * ( sss(:,:) - sb(:,:,1) ) / ( sn(:,:,1) + 1.e-20 )
712          ELSE                           ! NO Sea surface salinity damping
713             zsrp = 0.e0
714             erp(:,:) = 0.e0
715          ENDIF
716
717          ! freshwater fluxes terms
718          DO jj = 1, jpj
719             DO ji = 1, jpi
720                emp(ji, jj) =  emp(ji,jj) - zsumemp * tmask(ji,jj,1)
721                emps(ji, jj)=  emp(ji,jj) + erp(ji,jj)
722             END DO
723          END DO
724
725         IF( kt == nit000 .AND. lwp ) THEN
726            WRITE(numout,*)' ocesbc  : analytical formulation for gyre'
727            WRITE(numout,*)' ~~~~~~~ '
728            WRITE(numout,*)'           nyear      = ', nyear
729            WRITE(numout,*)'           ztime      = ',ztime
730            WRITE(numout,*)'           ztimemax1  = ',ztimemax1
731            WRITE(numout,*)'           ztimemin1  = ',ztimemin1
732            WRITE(numout,*)'           ztimemax2  = ',ztimemax2
733            WRITE(numout,*)'           ztimemin2  = ',ztimemin2
734            WRITE(numout,*)'           zyear0     = ', zyear0
735            WRITE(numout,*)'           zmonth0    = ', zmonth0
736            WRITE(numout,*)'           zday0      = ', zday0
737            WRITE(numout,*)'           zday_year0 = ',zday_year0
738            WRITE(numout,*)'           raajj      = ', raajj
739            WRITE(numout,*)'           zemp_S     = ',zemp_S 
740            WRITE(numout,*)'           zemp_N     = ',zemp_N 
741            WRITE(numout,*)'           zemp_sais  = ',zemp_sais
742            WRITE(numout,*)'           zTstar     = ',zTstar 
743            WRITE(numout,*)'           zsumemp    = ',zsumemp
744            WRITE(numout,*)'           zsurf      = ',zsurf 
745            WRITE(numout,*)'           dqdt0      = ',dqdt0
746            WRITE(numout,*)'           deds0      = ',deds0
747            WRITE(numout,*)'           zconv      = ',zconv 
748         ENDIF
749
750      ELSE
751
752         ! Constant surface fluxes
753
754         IF( kt == nit000 ) THEN
755            IF(lwp) THEN
756               WRITE(numout,*)' '
757               WRITE(numout,*)' ocesbc  : Constant surface fluxes read in namelist'
758               WRITE(numout,*)' ~~~~~~~ '
759               WRITE(numout,*)'           Namelist namflx: set the constant flux values'
760               WRITE(numout,*)'              net heat flux          q0   = ', q0  , ' W/m2'
761               WRITE(numout,*)'              solar heat flux        qsr0 = ', qsr0, ' W/m2'
762               WRITE(numout,*)'              net heat flux          emp0 = ', emp0, ' W/m2'
763            ENDIF
764
765            qt    (:,:) = q0
766            qsr   (:,:) = qsr0
767            emp   (:,:) = emp0
768            emps  (:,:) = emp0
769            qrp   (:,:) = 0.e0
770            erp   (:,:) = 0.e0
771
772            runoff(:,:) = 0.e0
773         ENDIF
774
775      ENDIF
776
777   END SUBROUTINE oce_sbc
778
779# endif
780#endif
781
782#if defined key_dtasal
783   !!----------------------------------------------------------------------
784   !!   'key_dtasal'                                          salinity data
785   !!----------------------------------------------------------------------
786   SUBROUTINE oce_sbc_dmp
787      !!---------------------------------------------------------------------
788      !!                   ***  ROUTINE oce_sbc_dmp  ***
789      !!                   
790      !! ** Purpose : Computation of internal and evaporation damping terms
791      !!        for ocean surface boundary conditions
792      !!
793      !! History :
794      !!   9.0  !  04-01  (G. Madec, C. Ethe)  Original code
795      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
796      !!----------------------------------------------------------------------
797      !! * Local declarations
798      INTEGER ::   ji, jj                   ! dummy loop indices
799      REAL(wp), DIMENSION(jpi,jpj)  :: zsss, zfreeze
800      REAL(wp) ::   zerp, zsrp
801!CT      CHARACTER (len=71) :: charout
802!CT#if ! defined key_dynspg_rl
803!CT      REAL(wp) ::   zwei
804!CT      REAL(wp) ::   zerpplus(jpi,jpj), zerpminus(jpi,jpj)
805!CT      REAL(wp) ::   zplus, zminus, zadefi
806!CT# if defined key_tradmp
807!CT      INTEGER jk
808!CT      REAL(wp), DIMENSION(jpi,jpj) ::   zstrdmp
809!CT# endif
810!CT#endif
811      !!----------------------------------------------------------------------
812
813#if defined key_ice_lim
814      ! sea ice indicator (1 or 0)
815      DO jj = 1, jpj
816         DO ji = 1, jpi
817            freezn(ji,jj) = MAX(0., SIGN(1., freeze(ji,jj)-rsmall) )
818         END DO
819      END DO
820      zsss   (:,:) = sss_io(:,:)
821      zfreeze(:,:) = freezn(:,:)
822#else
823      zsss   (:,:) = sb    (:,:,1)
824      zfreeze(:,:) = freeze(:,:)
825#endif
826
827      ! Initialisation
828      ! --------------
829      ! Restoring coefficients on SST and SSS   
830      zsrp = dqdt0 * ro0cpr * rauw   ! (Kg/m2/s)
831
832#if ! defined key_dynspg_rl
833      ! Free-surface
834         
835      ! Internal damping
836!CT# if defined key_tradmp
837!CT      ! Vertical mean of dampind trend (computed in tradmp module)
838!CT      zstrdmp(:,:) = 0.e0
839!CT      DO jk = 1, jpk
840!CT         zstrdmp(:,:) = zstrdmp(:,:) + strdmp(:,:,jk) * fse3t(:,:,jk)
841!CT      END DO
842!CT      ! volume flux associated to internal damping to climatology
843!CT      dmp(:,:) = zstrdmp(:,:) * rauw / ( zsss(:,:) + 1.e-20 )
844!CT# else
845      dmp(:,:) = 0.e0            ! No internal damping
846!CT# endif
847     
848      !   evaporation damping term ( Surface restoring )
849!CT      zerpplus (:,:) = 0.e0
850!CT      zerpminus(:,:) = 0.e0
851!CT      zplus  =  15. / rday
852!CT      zminus = -15. / rday
853     
854      DO jj = 1, jpj
855         DO ji = 1, jpi
856            zerp = ( 1. - 2.*upsrnfh(ji,jj) ) * zsrp   &
857               & * ( zsss(ji,jj) - s_dta(ji,jj,1) )     &
858               & / ( zsss(ji,jj) + 1.e-20        )
859           
860!CT            zerp = MIN( zerp, zplus  )
861!CT            zerp = MAX( zerp, zminus )
862            erp(ji,jj) = zerp
863!CT            zerpplus (ji,jj) = MAX( erp(ji,jj), 0.e0 )
864!CT            zerpminus(ji,jj) = MIN( erp(ji,jj), 0.e0 )
865         END DO
866      END DO
867
868!CT      aplus  = 0.e0
869!CT      aminus = 0.e0
870!CT     
871!CT      IF( nbit_cmp == 1) THEN
872!CT         
873!CT         IF(ln_ctl)   THEN
874!CT            WRITE(charout,FMT="('oce_sbc_dmp : a+ = ',D23.16, ' a- = ',D23.16)") aplus, aminus
875!CT            CALL prt_ctl_info(charout)
876!CT         ENDIF
877!CT         erp(:,:) = 0.e0
878!CT
879!CT      ELSE
880!CT
881!CT         DO jj = 1, jpj
882!CT            DO ji = 1, jpi
883!CT               zwei   = e1t(ji,jj) * e2t(ji,jj) * tmask_i(ji,jj)
884!CT               aplus  = aplus  + zerpplus (ji,jj) * zwei
885!CT               aminus = aminus - zerpminus(ji,jj) * zwei
886!CT            END DO
887!CT         END DO
888!CT         IF( lk_mpp )   CALL mpp_sum( aplus  )   ! sums over the global domain
889!CT         IF( lk_mpp )   CALL mpp_sum( aminus )
890!CT
891!CT         IF(ln_ctl)   THEN
892!CT            WRITE(charout,FMT="('oce_sbc_dmp : a+ = ',D23.16, ' a- = ',D23.16)") aplus, aminus
893!CT            CALL prt_ctl_info(charout)
894!CT         ENDIF
895!CT
896!CT         zadefi = MIN( aplus, aminus )
897!CT         IF( zadefi == 0.e0 ) THEN
898!CT            erp(:,:) = 0.e0
899!CT         ELSE
900!CT            erp(:,:) = zadefi * ( zerpplus(:,:) / aplus + zerpminus(:,:) / aminus )
901!CT         ENDIF
902!CT
903!CT      END IF
904#else
905      ! Rigid-lid (emp=emps=E-P-R+Erp)
906
907      erp(:,:) = ( 1. - zfreeze(:,:) ) * zsrp    &   ! surface restoring term
908         &     * ( zsss(:,:) - s_dta(:,:,1) )     &
909         &     / ( zsss(:,:) + 1.e-20      )
910#endif
911
912   END SUBROUTINE oce_sbc_dmp
913
914#else
915   !!----------------------------------------------------------------------
916   !!   Dummy routine                                      NO salinity data
917   !!----------------------------------------------------------------------
918   SUBROUTINE oce_sbc_dmp         ! Dummy routine
919      WRITE(*,*) 'oce_sbc_dmp: you should not have seen that print! error?'
920   END SUBROUTINE oce_sbc_dmp
921#endif
922
923   !!======================================================================
924END MODULE ocesbc
Note: See TracBrowser for help on using the repository browser.