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

Last change on this file since 556 was 531, checked in by opalod, 18 years ago

nemo_v1_update_75 : CT : enables bit comparison between single and multiple processor runs adding nbit_cmp namelist parameter

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