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

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

CL : Add CVS Header and CeCILL licence information

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