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

Last change on this file since 156 was 153, checked in by opalod, 20 years ago

CL + CT: BUGFIX098: Add missing declaration for kt and initialization of ztrp in coupling mode (key_coupled)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 33.2 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         END DO
373      END DO
374
375   END SUBROUTINE oce_sbc
376
377# elif defined key_flx_bulk_monthly || defined key_flx_bulk_daily || defined key_flx_forced_daily
378      !!-------------------------------------------------------------------------
379      !!   'key_flx_bulk_monthly' or 'key_flx_bulk_daily' or        bulk formulea
380      !!   'key_flx_forced_daily'                                or no bulk case
381      !!-------------------------------------------------------------------------
382
383   SUBROUTINE oce_sbc( kt )
384      !!---------------------------------------------------------------------
385      !!                   ***  ROUTINE oce_sbc  ***
386      !!                   
387      !! ** Purpose :   Ocean surface boundary conditions in forced mode
388      !!      using either flux or bulk formulation.
389      !!
390      !! History :
391      !!   1.0  !  99-11  (M. Imbard)  Original code
392      !!        !  01-03  (D. Ludicone, E. Durand, G. Madec) free surf.
393      !!   2.0  !  02-09  (G. Madec, C. Ethe)  F90: Free form and module
394      !!----------------------------------------------------------------------
395      !! * Modules used
396      USE daymod
397#if ! defined key_dtasst
398      USE dtasst, ONLY : rclice
399#endif
400#if defined key_flx_bulk_monthly || defined key_flx_bulk_daily
401      USE blk_oce
402#endif
403#if defined key_flx_forced_daily
404      USE flx_oce
405#endif
406
407      !! * arguments
408      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
409
410      !! * local declarations
411      INTEGER ::   ji, jj        ! dummy loop arguments
412      INTEGER ::   i15, ifreq             !     
413      REAL(wp) ::  zxy
414      REAL(wp) ::  zsice, zqri, zqrp, ztdta, zqrj
415      REAL(wp) ::  zq, zqi, zhemis, ztrp
416      REAL(wp), DIMENSION(jpi,jpj) :: zeri, zerps, ziclim
417      REAL(wp), DIMENSION(jpi,jpj) :: zqt, zqsr, zemp 
418      !!----------------------------------------------------------------------
419 
420      ! 1. initialization to zero at kt = nit000
421      ! ---------------------------------------
422     
423      IF( kt == nit000 ) THEN     
424         qsr    (:,:) = 0.e0
425         freeze (:,:) = 0.e0
426         qt     (:,:) = 0.e0
427         q      (:,:) = 0.e0
428         qrp    (:,:) = 0.e0
429         emp    (:,:) = 0.e0
430         emps   (:,:) = 0.e0
431         erp    (:,:) = 0.e0
432#if defined key_dynspg_fsc 
433         dmp    (:,:) = 0.e0
434#endif
435      ENDIF
436
437#if defined key_flx_bulk_monthly || defined key_flx_bulk_daily
438      ifreq      = nfbulk
439      zqt (:,:)  = qsr_oce(:,:) + qnsr_oce(:,:)
440      zqsr(:,:)  = qsr_oce(:,:)
441      zemp(:,:)  = evap(:,:) - tprecip(:,:)
442#endif 
443
444#if defined key_flx_forced_daily
445      ifreq      = 1
446      zqt (:,:)  = p_qt (:,:)
447      zqsr(:,:)  = p_qsr(:,:)
448      zemp(:,:)  = p_emp(:,:)
449#endif
450
451      IF( MOD( kt-1, ifreq) == 0 ) THEN
452         ! Computation of internal and evaporation damping terms       
453         CALL oce_sbc_dmp
454
455         ztrp = -40.             ! restoring terme for temperature (w/m2/k)   
456         zsice = - 0.04 / 0.8    ! ratio of isohaline compressibility over isotherme compressibility
457                                 ! ( d rho / dt ) / ( d rho / ds )      ( s = 34, t = -1.8 )
458         ! Flux computation
459         DO jj = 1, jpj
460            DO ji = 1, jpi     
461               ! climatological ice
462#if defined key_dtasst 
463               ziclim(ji,jj) = FLOAT( NINT( rclice(ji,jj,1) ) )
464#else
465               ! tested only with key key_dtasst (A. Lazar 07/2001)
466               ! this loop in CASE of interpolation of monthly rclice
467               i15           = INT( 2.* FLOAT(nday) / (FLOAT( nobis(nmonth) ) + 0.5) )
468               zxy           = FLOAT(nday) / FLOAT(nobis(nmonth)) + 0.5 - i15
469               ziclim(ji,jj) = FLOAT( NINT( (1-zxy) * rclice(ji,jj,1) + zxy  * rclice(ji,jj,2) ) )
470#endif
471
472               ! avoid surfreezing point           
473               tn(ji,jj,1) = MAX( tn(ji,jj,1), fzptn(ji,jj) )
474
475               ! hemisphere indicator (=1 north, =-1 south)           
476               zhemis = FLOAT( isign(1, mjg(jj) - (jpjdta/2+1) ) )
477
478               ! restoring temperature (ztdta >= to local freezing temperature)           
479#if defined key_dtasst
480               ztdta = MAX( sst(ji,jj),    fzptn(ji,jj) )
481#else
482               ztdta = MAX( t_dta(ji,jj,1), fzptn(ji,jj) )
483#endif
484
485               ! a) net downward radiative flux qsr()           
486               qsr(ji,jj) = zqsr(ji,jj) * tmask(ji,jj,1)
487
488               ! b) heat flux damping term qrp()
489               ! - gamma*(t-tlevitus) if no  climatological ice (ziclim=0)
490               ! - gamma*(t-(tgel-1.))  if climatological ice and no opa ice   (ziclim=1 zicopa=0)
491               ! - gamma*min(0,t-tgel) if climatological and opa ice (ziclim=1 zicopa=1)
492
493               zqrp = ztrp * ( tb(ji,jj,1) - ztdta )
494               zqri = ztrp * ( tb(ji,jj,1) - ( fzptn(ji,jj) - 1.) )
495               zqrj = ztrp * MIN( 0., tb(ji,jj,1) - fzptn(ji,jj) )
496               qrp(ji,jj) = ( (1. - ziclim(ji,jj)) * zqrp   &
497                  + ziclim(ji,jj)  * ( ( 1 - freeze(ji,jj)) * zqri   &
498                  + freeze(ji,jj)  * zqrj ) ) * tmask(ji,jj,1)
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          ztimemin = ztimemax + rjjhh * raajj / 2         ! rjjhh * raajj / 4 = 1 seasonal cycle in hours
654          ! amplitudes
655          zemp_S   = 0.7       ! intensity of COS in the South
656          zemp_N   = 0.8       ! intensity of COS in the North
657          zemp_sais= 0.1
658          zTstar   = 28.3      ! intemsity from 28.3 à -5°
659          zcos_sais = COS( (ztime - ztimemax) / (ztimemin - ztimemax) * rpi )    ! 1/2 period between 21th June and 21th December
660          ztrp= - 40.          ! retroaction term (W/m2/K)
661          zconv = 3.16e-5      ! convert 1m/yr->3.16e-5mm/s
662          DO jj = 1, jpj
663            DO ji = 1, jpi
664              t_star (ji,jj) = zTstar * ( 1 + 1. / 50. * zcos_sais )   &                  ! domain from 15° to 50°
665                 &             * COS( rpi * (gphit(ji,jj) - 5.)        &                  ! between 27 and 28 °C at 15N, -3
666                 &                        / (53.5 * ( 1 + 11 / 53.5 * zcos_sais ) * 2.) ) ! and 13°C at 50N 53.5 + or - 11
667                 !                                                                        ! = 1/4 period :
668                 !                                                                        ! 64.5 in summer, 42.5 in winter
669              qt  (ji,jj) = ztrp * ( tb(ji,jj,1) - t_star(ji,jj) )
670              IF( gphit(ji,jj) >= 14.845 .AND. 37.2 >= gphit(ji,jj)) THEN
671                 emp  (ji,jj) =   zemp_S * zconv   &
672                    &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (24.6 - 37.2) )  & ! zero at 37.8°, max at 24.6°
673                    &         * ( 1 - zemp_sais / zemp_S * zcos_sais)
674                 emps (ji,jj) = emp  (ji,jj)
675              ELSE
676                 emp (ji,jj) =  - zemp_N * zconv   &
677                    &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (46.8 - 37.2) )  & ! zero at 37.8°, max at 46.8°
678                    &         * ( 1 - zemp_sais / zemp_N * zcos_sais )
679                  emps (ji,jj) = emp  (ji,jj)
680              ENDIF
681              qsr (ji,jj) =  230 * COS( 3.1415 * ( gphit(ji,jj) - 23.5 * zcos_sais ) / ( 0.9 * 180 ) )   ! 23.5° : tropics
682            END DO
683          END DO
684          ! compute the emp flux such as its integration on the whole domain and at each time be zero
685          zsumemp = 0.
686          zsurf = 0.
687          DO jj = 1, jpj
688            DO ji = 1, jpi
689              zsumemp = zsumemp + emp(ji, jj) * tmask(ji, jj,  1)
690              zsurf =  zsurf +  tmask(ji, jj,  1)
691            END DO
692          END DO
693
694          IF( lk_mpp )   CALL mpp_sum( zsumemp  )       ! sum over the global domain
695          IF( lk_mpp )   CALL mpp_sum( zsurf    )       ! sum over the global domain
696
697          IF( nbench /= 0 ) THEN
698             ! Benchmark GYRE configuration (to allow the bit to bit comparison between Mpp/Mono case)
699             zsumemp = 0.e0
700          ELSE
701             ! Default GYRE configuration
702             zsumemp = zsumemp / zsurf
703          ENDIF
704          DO jj = 1, jpj
705            DO ji = 1, jpi
706              emp(ji, jj)=  emp(ji, jj) - zsumemp * tmask(ji, jj,  1)
707            END DO
708          END DO
709
710      ELSE
711
712         IF( kt == nit000 ) THEN
713
714            ! Read Namelist namflx : surface thermohaline fluxes
715            ! --------------------
716            REWIND ( numnam )
717            READ   ( numnam, namflx )
718   
719            IF(lwp) WRITE(numout,*)' '
720            IF(lwp) WRITE(numout,*)' ocesbc  : Constant surface fluxes read in namelist'
721            IF(lwp) WRITE(numout,*)' ~~~~~~~ '
722            IF(lwp) WRITE(numout,*)'           Namelist namflx: set the constant flux values'
723            IF(lwp) WRITE(numout,*)'              net heat flux          q0   = ', q0  , ' W/m2'
724            IF(lwp) WRITE(numout,*)'              solar heat flux        qsr0 = ', qsr0, ' W/m2'
725            IF(lwp) WRITE(numout,*)'              net heat flux          emp0 = ', emp0, ' W/m2'
726
727            qt    (:,:) = q0
728            qsr   (:,:) = qsr0
729            q     (:,:) = q0
730            emp   (:,:) = emp0
731            emps  (:,:) = emp0
732            qrp   (:,:) = 0.e0
733            erp   (:,:) = 0.e0
734   
735            runoff(:,:) = 0.e0
736         ENDIF
737      ENDIF
738
739   END SUBROUTINE oce_sbc
740
741# endif
742#endif
743
744#if defined key_dtasal
745   !!----------------------------------------------------------------------
746   !!   'key_dtasal'                                          salinity data
747   !!----------------------------------------------------------------------
748   SUBROUTINE oce_sbc_dmp
749      !!---------------------------------------------------------------------
750      !!                   ***  ROUTINE oce_sbc_dmp  ***
751      !!                   
752      !! ** Purpose : Computation of internal and evaporation damping terms
753      !!        for ocean surface boundary conditions
754      !!
755      !! History :
756      !!   9.0  !  04-01  (G. Madec, C. Ethe)  Original code
757      !!----------------------------------------------------------------------
758      !! * Local declarations
759      INTEGER ::   ji, jj                   ! dummy loop indices
760      REAL(wp), DIMENSION(jpi,jpj)  :: zsss, zfreeze
761      REAL(wp) ::   zerp, ztrp, zsrp
762#if defined key_dynspg_fsc
763      REAL(wp) ::   zwei
764      REAL(wp) ::   zerpplus(jpi,jpj), zerpminus(jpi,jpj)
765      REAL(wp) ::   zplus, zminus, zadefi
766# if defined key_tradmp
767      INTEGER jk
768      REAL(wp), DIMENSION(jpi,jpj) ::   zstrdmp
769# endif
770#endif
771      !!----------------------------------------------------------------------
772
773#if defined key_ice_lim
774      ! sea ice indicator (1 or 0)
775      DO jj = 1, jpj
776         DO ji = 1, jpi
777            freezn(ji,jj) = MAX(0., SIGN(1., freeze(ji,jj)-rsmall) )
778         END DO
779      END DO
780      zsss   (:,:) = sss_io(:,:)
781      zfreeze(:,:) = freezn(:,:)
782#else
783      zsss   (:,:) = sn    (:,:,1)
784      zfreeze(:,:) = freeze(:,:)
785#endif
786
787      ! Initialisation
788      ! --------------
789      ! Restoring coefficients on SST and SSS   
790      IF( lk_cpl ) THEN
791         ztrp = 0.e0
792         zsrp = 0.e0
793      ELSE
794         ztrp  = -40.                   ! (W/m2/K)   
795         zsrp  = ztrp * ro0cpr * rauw   ! (Kg/m2/s2)
796      ENDIF
797
798#if defined key_dynspg_fsc 
799      ! Free-surface
800         
801      ! Internal damping
802# if defined key_tradmp
803      ! Vertical mean of dampind trend (computed in tradmp module)
804      zstrdmp(:,:) = 0.e0
805      DO jk = 1, jpk
806         zstrdmp(:,:) = zstrdmp(:,:) + strdmp(:,:,jk) * fse3t(:,:,jk)
807      END DO
808      ! volume flux associated to internal damping to climatology
809!!ibu dmp(:,:) = zstrdmp(:,:) * rauw / ( zsss(:,:) + rsmall )
810      dmp(:,:) = zstrdmp(:,:) * rauw / ( zsss(:,:) + 1.e-20 )
811# else
812      dmp(:,:) = 0.e0            ! No internal damping
813# endif
814     
815      !   evaporation damping term ( Surface restoring )
816      zerpplus (:,:) = 0.e0
817      zerpminus(:,:) = 0.e0
818      zplus  =  15. / rday
819      zminus = -15. / rday
820     
821      DO jj = 1, jpj
822         DO ji = 1, jpi
823            zerp = ( 1. - 2.*upsrnfh(ji,jj) ) * zsrp   &
824               & * ( zsss(ji,jj) - s_dta(ji,jj,1) )     &
825               & / ( zsss(ji,jj) + 1.e-20        )
826!ib            & / ( zsss(ji,jj) + rsmall        )
827           
828            zerp = MIN( zerp, zplus  )
829            zerp = MAX( zerp, zminus )
830            erp(ji,jj) = zerp
831            zerpplus (ji,jj) = MAX( erp(ji,jj), 0.e0 )
832            zerpminus(ji,jj) = MIN( erp(ji,jj), 0.e0 )
833         END DO
834      END DO
835
836      aplus  = 0.e0
837      aminus = 0.e0
838      DO jj = 1, jpj
839         DO ji = 1, jpi
840            zwei   = e1t(ji,jj) * e2t(ji,jj) * tmask_i(ji,jj)
841            aplus  = aplus  + zerpplus (ji,jj) * zwei
842            aminus = aminus - zerpminus(ji,jj) * zwei
843         END DO
844      END DO
845      IF( lk_mpp )   CALL mpp_sum( aplus  )   ! sums over the global domain
846      IF( lk_mpp )   CALL mpp_sum( aminus )
847      IF(l_ctl)   WRITE(numout,*) ' oce_sbc_dmp : a+ = ', aplus, ' a- = ', aminus
848
849      zadefi = MIN( aplus, aminus )
850      IF( zadefi == 0.e0 ) THEN
851         erp(:,:) = 0.e0
852      ELSE
853         erp(:,:) = zadefi * ( zerpplus(:,:) / aplus + zerpminus(:,:) / aminus )
854      ENDIF
855#else
856      ! Rigid-lid (emp=emps=E-P-R+Erp)
857     
858      erp(:,:) = ( 1. - zfreeze(:,:) ) * zsrp    &   ! surface restoring term
859         &     * ( zsss(:,:) - s_dta(:,:,1) )     &
860         &     / ( zsss(:,:) + 1.e-20      )
861!ib      &     / ( zsss(:,:) + rsmall      )
862#endif
863
864   END SUBROUTINE oce_sbc_dmp
865
866#else
867   !!----------------------------------------------------------------------
868   !!   Dummy routine                                      NO salinity data
869   !!----------------------------------------------------------------------
870   SUBROUTINE oce_sbc_dmp         ! Dummy routine
871      WRITE(*,*) 'oce_sbc_dmp: you should not have seen that print! error?'
872   END SUBROUTINE oce_sbc_dmp
873#endif
874
875   !!======================================================================
876END MODULE ocesbc
Note: See TracBrowser for help on using the repository browser.