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

Last change on this file since 227 was 227, checked in by opalod, 19 years ago

CT : BUGFIX158 : - correction of inconsistencies when using bulk formulae without the sea-ice model

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