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 tags/nemo_dev_x3/NEMO/OPA_SRC/SBC – NEMO

source: tags/nemo_dev_x3/NEMO/OPA_SRC/SBC/ocesbc.F90 @ 105

Last change on this file since 105 was 105, checked in by cvs2svn, 20 years ago

This commit was manufactured by cvs2svn to create tag 'nemo_dev_x3'.

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