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

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

nemo_v1_bugfix_003 : CT : when not using sea-ice the damping term qrp(:,:) must be computed only if not using the monthly bulk neither the daily bulk

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