source: trunk/NEMO/OPA_SRC/istate.F90 @ 900

Last change on this file since 900 was 900, checked in by rblod, 13 years ago

Update 1D configuration according to SBC and LIM3, see ticket #117

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 26.9 KB
Line 
1MODULE istate
2   !!======================================================================
3   !!                     ***  MODULE  istate  ***
4   !! Ocean state   :  initial state setting
5   !!=====================================================================
6   !! History :   4.0  !  89-12  (P. Andrich)  Original code
7   !!             5.0  !  91-11  (G. Madec)  rewritting
8   !!             6.0  !  96-01  (G. Madec)  terrain following coordinates
9   !!             8.0  !  01-09  (M. Levy, M. Ben Jelloul)  istate_eel
10   !!             8.0  !  01-09  (M. Levy, M. Ben Jelloul)  istate_uvg
11   !!             9.0  !  03-08  (G. Madec)  F90: Free form, modules
12   !!             9.0  !  03-09  (G. Madec, C. Talandier)  add EEL R5
13   !!             9.0  !  04-05  (A. Koch-Larrouy)  istate_gyre
14   !!             9.0  !  06-07  (S. Masson)  distributed restart using iom
15   !!----------------------------------------------------------------------
16
17   !!----------------------------------------------------------------------
18   !!   istate_init   : initial state setting
19   !!   istate_tem    : analytical profile for initial Temperature
20   !!   istate_sal    : analytical profile for initial Salinity
21   !!   istate_eel    : initial state setting of EEL R5 configuration
22   !!   istate_gyre   : initial state setting of GYRE configuration
23   !!   istate_uvg    : initial velocity in geostropic balance
24   !!----------------------------------------------------------------------
25   USE oce             ! ocean dynamics and active tracers
26   USE dom_oce         ! ocean space and time domain
27   USE daymod          !
28   USE ldftra_oce      ! ocean active tracers: lateral physics
29   USE zdf_oce         ! ocean vertical physics
30   USE phycst          ! physical constants
31   USE wzvmod          ! verctical velocity               (wzv     routine)
32   USE dtatem          ! temperature data                 (dta_tem routine)
33   USE dtasal          ! salinity data                    (dta_sal routine)
34   USE restart         ! ocean restart                   (rst_read routine)
35   USE solisl          ! ???
36   USE in_out_manager  ! I/O manager
37   USE iom
38   USE c1d             ! re-initialization of u-v mask for the 1D configuration
39   USE zpshde          ! partial step: hor. derivative (zps_hde routine)
40   USE eosbn2          ! equation of state            (eos bn2 routine)
41   USE domvvl          ! varying vertical mesh
42   USE dynspg_oce      ! pressure gradient schemes
43   USE dynspg_flt      ! pressure gradient schemes
44   USE dynspg_exp      ! pressure gradient schemes
45   USE dynspg_ts       ! pressure gradient schemes
46   
47   IMPLICIT NONE
48   PRIVATE
49
50   PUBLIC   istate_init   ! routine called by step.F90
51
52   !! * Substitutions
53#  include "domzgr_substitute.h90"
54#  include "vectopt_loop_substitute.h90"
55   !!----------------------------------------------------------------------
56   !!   OPA 9.0 , LOCEAN-IPSL (2006)
57   !! $Id$
58   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
59   !!----------------------------------------------------------------------
60
61CONTAINS
62
63   SUBROUTINE istate_init
64      !!----------------------------------------------------------------------
65      !!                   ***  ROUTINE istate_init  ***
66      !!
67      !! ** Purpose :   Initialization of the dynamics and tracer fields.
68      !!----------------------------------------------------------------------
69      USE eosbn2          ! eq. of state, Brunt Vaisala frequency (eos     routine)
70
71      IF(lwp) WRITE(numout,*)
72      IF(lwp) WRITE(numout,*) 'istate_ini : Initialization of the dynamics and tracers'
73      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
74
75      rhd  (:,:,:) = 0.e0
76      rhop (:,:,:) = 0.e0
77      rn2  (:,:,:) = 0.e0 
78
79      IF( ln_rstart ) THEN                    ! Restart from a file
80         !                                    ! -------------------
81         neuler = 1                              ! Set time-step indicator at nit000 (leap-frog)
82         CALL rst_read                           ! Read the restart file
83         CALL day( -1 )                          ! model calendar (using both namelist and restart infos)
84      ELSE
85         !                                    ! Start from rest
86         !                                    ! ---------------
87         neuler = 0                              ! Set time-step indicator at nit000 (euler forward)
88         adatrj = 0._wp
89         CALL day( -1 )                          ! model calendar (using namelist infos)
90         numror = 0                              ! define numror = 0 -> no restart file to read
91         !                                       ! Initialization of ocean to zero
92         !     before fields       !       now fields         
93             ub   (:,:,:) = 0.e0   ;   un   (:,:,:) = 0.e0   ; sshb(:,:) = 0.e0
94             vb   (:,:,:) = 0.e0   ;   vn   (:,:,:) = 0.e0   ; sshn(:,:) = 0.e0
95             rotb (:,:,:) = 0.e0   ;   rotn (:,:,:) = 0.e0
96             hdivb(:,:,:) = 0.e0   ;   hdivn(:,:,:) = 0.e0
97         !
98         IF( cp_cfg == 'eel' ) THEN
99            CALL istate_eel                      ! EEL   configuration : start from pre-defined
100            !                                    !                       velocity and thermohaline fields
101         ELSEIF( cp_cfg == 'gyre' ) THEN         
102            CALL istate_gyre                     ! GYRE  configuration : start from pre-defined temperature
103            !                                    !                       and salinity fields
104         ELSE
105            !                                    ! Other configurations: Initial temperature and salinity fields
106#if defined key_dtatem
107            CALL dta_tem( nit000 )                  ! read 3D temperature data
108            tb(:,:,:) = t_dta(:,:,:)                ! use temperature data read
109            tn(:,:,:) = t_dta(:,:,:)
110#else
111            IF(lwp) WRITE(numout,*)                 ! analytical temperature profile
112            IF(lwp) WRITE(numout,*)'             Temperature initialization using an analytic profile'
113            CALL istate_tem
114#endif
115#if defined key_dtasal
116            CALL dta_sal( nit000 )                  ! read 3D salinity data
117            sb(:,:,:) = s_dta(:,:,:)                ! use salinity data read
118            sn(:,:,:) = s_dta(:,:,:)
119#else
120            ! No salinity data
121            IF(lwp)WRITE(numout,*)                  ! analytical salinity profile
122            IF(lwp)WRITE(numout,*)'             Salinity initialisation using a constant value'
123            CALL istate_sal
124#endif
125         ENDIF
126
127         CALL eos( tb, sb, rhd, rhop )        ! before potential and in situ densities
128         
129         IF( ln_zps .AND. .NOT. lk_c1d )   &
130            &             CALL zps_hde( nit000, tb, sb, rhd,  &  ! Partial steps: before Horizontal DErivative
131            &                                  gtu, gsu, gru, &  ! of t, s, rd at the bottom ocean level
132            &                                  gtv, gsv, grv )
133         
134      ENDIF
135
136      IF( lk_vvl ) THEN
137         ! read free surface arrays in restart file
138         IF( ln_rstart ) THEN
139            IF( lk_dynspg_flt )   CALL flt_rst( nit000, 'READ' )      ! read or initialize the following fields
140            !                                                         ! gcx, gcxb, sshb, sshn
141            IF( lk_dynspg_ts  )   CALL ts_rst ( nit000, 'READ' )      ! read or initialize the following fields
142            !                                                         ! sshb, sshn, sshb_b, sshn_b, un_b, vn_b
143            IF( lk_dynspg_exp )   CALL exp_rst( nit000, 'READ' )      ! read or initialize the following fields
144            !                                                         ! sshb, sshn
145         ENDIF
146         !
147         IF( .NOT. lk_dynspg_flt ) sshbb(:,:) = sshb(:,:)
148         !
149         CALL dom_vvl               ! ssh init and vertical grid update
150
151         CALL eos_init              ! usefull to get the equation state type neos parameter
152
153         CALL bn2( tb, sb, rn2 )    ! before Brunt Vaissala frequency
154
155         IF( .NOT. ln_rstart ) CALL wzv( nit000 ) 
156
157      ENDIF
158
159      !                                       ! Vertical velocity
160      !                                       ! -----------------
161
162      IF( .NOT. lk_vvl )    CALL wzv( nit000 )                         ! from horizontal divergence
163      !
164   END SUBROUTINE istate_init
165
166
167   SUBROUTINE istate_tem
168      !!---------------------------------------------------------------------
169      !!                  ***  ROUTINE istate_tem  ***
170      !!   
171      !! ** Purpose :   Intialization of the temperature field with an
172      !!      analytical profile or a file (i.e. in EEL configuration)
173      !!
174      !! ** Method  :   Use Philander analytic profile of temperature
175      !!
176      !! References :  Philander ???
177      !!----------------------------------------------------------------------
178      INTEGER :: ji, jj, jk
179      !!----------------------------------------------------------------------
180      !
181      IF(lwp) WRITE(numout,*)
182      IF(lwp) WRITE(numout,*) 'istate_tem : initial temperature profile'
183      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
184
185      DO jk = 1, jpk
186         DO jj = 1, jpj
187            DO ji = 1, jpi
188               tn(ji,jj,jk) = (  ( ( 7.5 - 0.*ABS(gphit(ji,jj))/30. )   &
189                  &               *( 1.-TANH((fsdept(ji,jj,jk)-80.)/30.) )   &
190                  &            + 10.*(5000.-fsdept(ji,jj,jk))/5000.)  ) * tmask(ji,jj,jk)
191               tb(ji,jj,jk) = tn(ji,jj,jk)
192          END DO
193        END DO
194      END DO
195
196      IF(lwp) CALL prizre( tn    , jpi   , jpj   , jpk   , jpj/2 ,   &
197         &                 1     , jpi   , 5     , 1     , jpk   ,   &
198         &                 1     , 1.    , numout                  )
199      !
200   END SUBROUTINE istate_tem
201
202
203   SUBROUTINE istate_sal
204      !!---------------------------------------------------------------------
205      !!                  ***  ROUTINE istate_sal  ***
206      !!
207      !! ** Purpose :   Intialize the salinity field with an analytic profile
208      !!
209      !! ** Method  :   Use to a constant value 35.5
210      !!             
211      !! ** Action  :   Initialize sn and sb
212      !!----------------------------------------------------------------------
213      REAL(wp) ::   zsal = 35.50_wp
214      !!----------------------------------------------------------------------
215
216      IF(lwp) WRITE(numout,*)
217      IF(lwp) WRITE(numout,*) 'istate_sal : initial salinity : ', zsal
218      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
219
220      sn(:,:,:) = zsal * tmask(:,:,:)
221      sb(:,:,:) = sn(:,:,:)
222     
223   END SUBROUTINE istate_sal
224
225
226   SUBROUTINE istate_eel
227      !!----------------------------------------------------------------------
228      !!                   ***  ROUTINE istate_eel  ***
229      !!
230      !! ** Purpose :   Initialization of the dynamics and tracers for EEL R5
231      !!      configuration (channel with or without a topographic bump)
232      !!
233      !! ** Method  : - set temprature field
234      !!              - set salinity field
235      !!              - set velocity field including horizontal divergence
236      !!                and relative vorticity fields
237      !!----------------------------------------------------------------------
238      USE eosbn2     ! eq. of state, Brunt Vaisala frequency (eos     routine)
239      USE divcur     ! hor. divergence & rel. vorticity      (div_cur routine)
240      USE iom
241 
242      INTEGER  ::   inum              ! temporary logical unit
243      INTEGER  ::   ji, jj, jk        ! dummy loop indices
244      INTEGER  ::   ijloc
245      REAL(wp) ::   zh1, zh2, zslope, zcst, zfcor   ! temporary scalars
246      REAL(wp) ::   zt1  = 15._wp,               &  ! surface temperature value (EEL R5)
247         &          zt2  =  5._wp,               &  ! bottom  temperature value (EEL R5)
248         &          zsal = 35.0_wp,              &  ! constant salinity (EEL R2, R5 and R6)
249         &          zueel = 0.1_wp                  ! constant uniform zonal velocity (EEL R5)
250# if ! defined key_dynspg_rl
251      REAL(wp), DIMENSION(jpiglo,jpjglo) ::   zssh  ! initial ssh over the global domain
252# endif
253      !!----------------------------------------------------------------------
254
255      SELECT CASE ( jp_cfg ) 
256         !                                              ! ====================
257         CASE ( 5 )                                     ! EEL R5 configuration
258            !                                           ! ====================
259
260            ! set temperature field with a linear profile
261            ! -------------------------------------------
262            IF(lwp) WRITE(numout,*)
263            IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: linear temperature profile'
264            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
265
266            zh1 = gdept_0(  1  )
267            zh2 = gdept_0(jpkm1)
268
269            zslope = ( zt1 - zt2 ) / ( zh1 - zh2 )
270            zcst   = ( zt1 * ( zh1 - zh2) - ( zt1 - zt2 ) * zh1 ) / ( zh1 - zh2 )
271
272            DO jk = 1, jpk
273               tn(:,:,jk) = ( zt2 + zt1 * exp( - fsdept(:,:,jk) / 1000 ) ) * tmask(:,:,jk)
274               tb(:,:,jk) = tn(:,:,jk)
275            END DO
276
277            IF(lwp) CALL prizre( tn    , jpi   , jpj   , jpk   , jpj/2 ,   &
278               &                 1     , jpi   , 5     , 1     , jpk   ,   &
279               &                 1     , 1.    , numout                  )
280
281            ! set salinity field to a constant value
282            ! --------------------------------------
283            IF(lwp) WRITE(numout,*)
284            IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: constant salinity field, S = ', zsal
285            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
286
287            sn(:,:,:) = zsal * tmask(:,:,:)
288            sb(:,:,:) = sn(:,:,:)
289     
290
291# if ! defined key_dynspg_rl
292            ! set the dynamics: U,V, hdiv, rot (and ssh if necessary)
293            ! ----------------
294            ! Start EEL5 configuration with barotropic geostrophic velocities
295            ! according the sshb and sshn SSH imposed.
296            ! we assume a uniform grid (hence the use of e1t(1,1) for delta_y)
297            ! we use the Coriolis frequency at mid-channel.   
298   
299            ub(:,:,:) = zueel * umask(:,:,:)
300            un(:,:,:) = ub(:,:,:)
301            ijloc = mj0(INT(jpjglo-1)/2)
302            zfcor = ff(1,ijloc)
303
304            DO jj = 1, jpjglo
305               zssh(:,jj) = - (FLOAT(jj)- FLOAT(jpjglo-1)/2.)*zueel*e1t(1,1)*zfcor/grav 
306            END DO
307
308            IF(lwp) THEN
309               WRITE(numout,*) ' Uniform zonal velocity for EEL R5:',zueel
310               WRITE(numout,*) ' Geostrophic SSH profile as a function of y:'
311               WRITE(numout,'(12(1x,f6.2))') zssh(1,:)
312            ENDIF
313
314            DO jj = 1, nlcj
315               DO ji = 1, nlci
316                  sshb(ji,jj) = zssh( mig(ji) , mjg(jj) ) * tmask(ji,jj,1)
317               END DO
318            END DO
319            sshb(nlci+1:jpi,      :   ) = 0.e0      ! set to zero extra mpp columns
320            sshb(      :   ,nlcj+1:jpj) = 0.e0      ! set to zero extra mpp rows
321
322            sshn(:,:) = sshb(:,:)                   ! set now ssh to the before value
323
324            IF( nn_rstssh /= 0 ) THEN 
325               nn_rstssh = 0                           ! hand-made initilization of ssh
326               CALL ctl_warn( 'istate_eel: force nn_rstssh = 0' )
327            ENDIF
328
329            ! horizontal divergence and relative vorticity (curl)
330            CALL div_cur( nit000 )
331
332            ! N.B. the vertical velocity will be computed from the horizontal divergence field
333            ! in istate by a call to wzv routine
334# endif
335
336
337            !                                     ! ==========================
338         CASE ( 2 , 6 )                           ! EEL R2 or R6 configuration
339            !                                     ! ==========================
340 
341            ! set temperature field with a NetCDF file
342            ! ----------------------------------------
343            IF(lwp) WRITE(numout,*)
344            IF(lwp) WRITE(numout,*) 'istate_eel : EEL R2 or R6: read initial temperature in a NetCDF file'
345            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
346
347            CALL iom_open ( 'eel.initemp', inum )
348            CALL iom_get ( inum, jpdom_data, 'initemp', tb ) ! read before temprature (tb)
349            CALL iom_close( inum )
350     
351            tn(:,:,:) = tb(:,:,:)                            ! set nox temperature to tb
352
353            IF(lwp) CALL prizre( tn    , jpi   , jpj   , jpk   , jpj/2 ,   &
354               &                 1     , jpi   , 5     , 1     , jpk   ,   &
355               &                 1     , 1.    , numout                  )
356
357
358            ! set salinity field to a constant value
359            ! --------------------------------------
360            IF(lwp) WRITE(numout,*)
361            IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: constant salinity field, S = ', zsal
362            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
363 
364            sn(:,:,:) = zsal * tmask(:,:,:)
365            sb(:,:,:) = sn(:,:,:)
366
367
368            IF( lk_isl ) THEN
369               ! Horizontal velocity : start from geostrophy (EEL config)
370               CALL eos( tn, sn, rhd )     ! now in situ density
371               CALL istate_uvg             ! compute geostrophic velocity
372
373               ! N.B. the vertical velocity will be computed from the horizontal divergence field
374               ! in istate by a call to wzv routine
375            ENDIF
376            !                                    ! ===========================
377         CASE DEFAULT                            ! NONE existing configuration
378            !                                    ! ===========================
379            WRITE(ctmp1,*) 'EEL with a ', jp_cfg,' km resolution is not coded'
380            CALL ctl_stop( ctmp1 )
381
382      END SELECT
383
384   END SUBROUTINE istate_eel
385
386
387   SUBROUTINE istate_gyre
388      !!----------------------------------------------------------------------
389      !!                   ***  ROUTINE istate_gyre  ***
390      !!
391      !! ** Purpose :   Initialization of the dynamics and tracers for GYRE
392      !!      configuration (double gyre with rotated domain)
393      !!
394      !! ** Method  : - set temprature field
395      !!              - set salinity field
396      !!----------------------------------------------------------------------
397      INTEGER :: ji, jj, jk  ! dummy loop indices
398      INTEGER            ::   inum          ! temporary logical unit
399      INTEGER, PARAMETER ::   ntsinit = 0   ! (0/1) (analytical/input data files) T&S initialization
400      !!----------------------------------------------------------------------
401
402      SELECT CASE ( ntsinit)
403
404      CASE ( 0 )                  ! analytical T/S profil deduced from LEVITUS
405         IF(lwp) WRITE(numout,*)
406         IF(lwp) WRITE(numout,*) 'istate_gyre : initial analytical T and S profil deduced from LEVITUS '
407         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
408
409         DO jk = 1, jpk
410            DO jj = 1, jpj
411               DO ji = 1, jpi
412                  tn(ji,jj,jk) = (  16. - 12. * TANH( (fsdept(ji,jj,jk) - 400) / 700 )         )   &
413                       &           * (-TANH( (500-fsdept(ji,jj,jk)) / 150 ) + 1) / 2               &
414                       &       + (      15. * ( 1. - TANH( (fsdept(ji,jj,jk)-50.) / 1500.) )       &
415                       &                - 1.4 * TANH((fsdept(ji,jj,jk)-100.) / 100.)               &   
416                       &                + 7.  * (1500. - fsdept(ji,jj,jk)) / 1500.             )   & 
417                       &           * (-TANH( (fsdept(ji,jj,jk) - 500) / 150) + 1) / 2
418                  tn(ji,jj,jk) = tn(ji,jj,jk) * tmask(ji,jj,jk)
419                  tb(ji,jj,jk) = tn(ji,jj,jk)
420
421                  sn(ji,jj,jk) =  (  36.25 - 1.13 * TANH( (fsdept(ji,jj,jk) - 305) / 460 )  )  &
422                     &              * (-TANH((500 - fsdept(ji,jj,jk)) / 150) + 1) / 2          &
423                     &          + (  35.55 + 1.25 * (5000. - fsdept(ji,jj,jk)) / 5000.         &
424                     &                - 1.62 * TANH( (fsdept(ji,jj,jk) - 60.  ) / 650. )       &
425                     &                + 0.2  * TANH( (fsdept(ji,jj,jk) - 35.  ) / 100. )       &
426                     &                + 0.2  * TANH( (fsdept(ji,jj,jk) - 1000.) / 5000.)    )  &
427                     &              * (-TANH((fsdept(ji,jj,jk) - 500) / 150) + 1) / 2 
428                  sn(ji,jj,jk) = sn(ji,jj,jk) * tmask(ji,jj,jk)
429                  sb(ji,jj,jk) = sn(ji,jj,jk)
430               END DO
431            END DO
432         END DO
433
434      CASE ( 1 )                  ! T/S data fields read in dta_tem.nc/data_sal.nc files
435         IF(lwp) WRITE(numout,*)
436         IF(lwp) WRITE(numout,*) 'istate_gyre : initial T and S read from dta_tem.nc/data_sal.nc files'
437         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
438         IF(lwp) WRITE(numout,*) '              NetCDF FORMAT'
439
440         ! Read temperature field
441         ! ----------------------
442         CALL iom_open ( 'data_tem', inum )
443         CALL iom_get ( inum, jpdom_data, 'votemper', tn ) 
444         CALL iom_close( inum )
445
446         tn(:,:,:) = tn(:,:,:) * tmask(:,:,:) 
447         tb(:,:,:) = tn(:,:,:)
448
449         ! Read salinity field
450         ! -------------------
451         CALL iom_open ( 'data_sal', inum )
452         CALL iom_get ( inum, jpdom_data, 'vosaline', sn ) 
453         CALL iom_close( inum )
454
455         sn(:,:,:)  = sn(:,:,:) * tmask(:,:,:) 
456         sb(:,:,:)  = sn(:,:,:)
457
458      END SELECT
459
460      IF(lwp) THEN
461         WRITE(numout,*)
462         WRITE(numout,*) '              Initial temperature and salinity profiles:'
463         WRITE(numout, "(9x,' level   gdept_0   temperature   salinity   ')" )
464         WRITE(numout, "(10x, i4, 3f10.2)" ) ( jk, gdept_0(jk), tn(2,2,jk), sn(2,2,jk), jk = 1, jpk )
465      ENDIF
466
467   END SUBROUTINE istate_gyre
468
469
470   SUBROUTINE istate_uvg
471      !!----------------------------------------------------------------------
472      !!                  ***  ROUTINE istate_uvg  ***
473      !!
474      !! ** Purpose :   Compute the geostrophic velocities from (tn,sn) fields
475      !!
476      !! ** Method  :   Using the hydrostatic hypothesis the now hydrostatic
477      !!      pressure is computed by integrating the in-situ density from the
478      !!      surface to the bottom.
479      !!                 p=integral [ rau*g dz ]
480      !!----------------------------------------------------------------------
481      USE dynspg          ! surface pressure gradient             (dyn_spg routine)
482      USE divcur          ! hor. divergence & rel. vorticity      (div_cur routine)
483      USE lbclnk          ! ocean lateral boundary condition (or mpp link)
484
485      INTEGER ::   ji, jj, jk        ! dummy loop indices
486      INTEGER ::   indic             ! ???
487      REAL(wp) ::   zmsv, zphv, zmsu, zphu, zalfg     ! temporary scalars
488      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zprn     ! workspace
489      !!----------------------------------------------------------------------
490
491      IF(lwp) WRITE(numout,*) 
492      IF(lwp) WRITE(numout,*) 'istate_uvg : Start from Geostrophy'
493      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
494
495      ! Compute the now hydrostatic pressure
496      ! ------------------------------------
497
498      zalfg = 0.5 * grav * rau0
499     
500      zprn(:,:,1) = zalfg * fse3w(:,:,1) * ( 1 + rhd(:,:,1) )       ! Surface value
501
502      DO jk = 2, jpkm1                                              ! Vertical integration from the surface
503         zprn(:,:,jk) = zprn(:,:,jk-1)   &
504            &         + zalfg * fse3w(:,:,jk) * ( 2. + rhd(:,:,jk) + rhd(:,:,jk-1) )
505      END DO 
506
507      ! Compute geostrophic balance
508      ! ---------------------------
509      DO jk = 1, jpkm1
510         DO jj = 2, jpjm1
511            DO ji = fs_2, fs_jpim1   ! vertor opt.
512               zmsv = 1. / MAX(  umask(ji-1,jj+1,jk) + umask(ji  ,jj+1,jk)   &
513                               + umask(ji-1,jj  ,jk) + umask(ji  ,jj  ,jk) , 1.  )
514               zphv = ( zprn(ji  ,jj+1,jk) - zprn(ji-1,jj+1,jk) ) * umask(ji-1,jj+1,jk) / e1u(ji-1,jj+1)   &
515                    + ( zprn(ji+1,jj+1,jk) - zprn(ji  ,jj+1,jk) ) * umask(ji  ,jj+1,jk) / e1u(ji  ,jj+1)   &
516                    + ( zprn(ji  ,jj  ,jk) - zprn(ji-1,jj  ,jk) ) * umask(ji-1,jj  ,jk) / e1u(ji-1,jj  )   &
517                    + ( zprn(ji+1,jj  ,jk) - zprn(ji  ,jj  ,jk) ) * umask(ji  ,jj  ,jk) / e1u(ji  ,jj  )
518               zphv = 1. / rau0 * zphv * zmsv * vmask(ji,jj,jk)
519
520               zmsu = 1. / MAX(  vmask(ji+1,jj  ,jk) + vmask(ji  ,jj  ,jk)   &
521                               + vmask(ji+1,jj-1,jk) + vmask(ji  ,jj-1,jk) , 1.  )
522               zphu = ( zprn(ji+1,jj+1,jk) - zprn(ji+1,jj  ,jk) ) * vmask(ji+1,jj  ,jk) / e2v(ji+1,jj  )   &
523                    + ( zprn(ji  ,jj+1,jk) - zprn(ji  ,jj  ,jk) ) * vmask(ji  ,jj  ,jk) / e2v(ji  ,jj  )   &
524                    + ( zprn(ji+1,jj  ,jk) - zprn(ji+1,jj-1,jk) ) * vmask(ji+1,jj-1,jk) / e2v(ji+1,jj-1)   &
525                    + ( zprn(ji  ,jj  ,jk) - zprn(ji  ,jj-1,jk) ) * vmask(ji  ,jj-1,jk) / e2v(ji  ,jj-1)
526               zphu = 1. / rau0 * zphu * zmsu * umask(ji,jj,jk)
527
528               ! Compute the geostrophic velocities
529               un(ji,jj,jk) = -2. * zphu / ( ff(ji,jj) + ff(ji  ,jj-1) )
530               vn(ji,jj,jk) =  2. * zphv / ( ff(ji,jj) + ff(ji-1,jj  ) )
531            END DO
532         END DO
533      END DO
534
535      IF(lwp) WRITE(numout,*) '         we force to zero bottom velocity'
536
537      ! Susbtract the bottom velocity (level jpk-1 for flat bottom case)
538      ! to have a zero bottom velocity
539
540      DO jk = 1, jpkm1
541         un(:,:,jk) = ( un(:,:,jk) - un(:,:,jpkm1) ) * umask(:,:,jk)
542         vn(:,:,jk) = ( vn(:,:,jk) - vn(:,:,jpkm1) ) * vmask(:,:,jk)
543      END DO
544
545      CALL lbc_lnk( un, 'U', -1. )
546      CALL lbc_lnk( vn, 'V', -1. )
547     
548      ub(:,:,:) = un(:,:,:)
549      vb(:,:,:) = vn(:,:,:)
550     
551      ! WARNING !!!!!
552      ! after initializing u and v, we need to calculate the initial streamfunction bsf.
553      ! Otherwise, only the trend will be computed and the model will blow up (inconsistency).
554      ! to do that, we call dyn_spg with a special trick:
555      ! we fill ua and va with the velocities divided by dt, and the streamfunction will be brought to the
556      ! right value assuming the velocities have been set up in one time step.
557      ! we then set bsfd to zero (first guess for next step is d(psi)/dt = 0.)
558      !  sets up s false trend to calculate the barotropic streamfunction.
559
560      ua(:,:,:) = ub(:,:,:) / rdt
561      va(:,:,:) = vb(:,:,:) / rdt
562
563      ! calls dyn_spg. we assume euler time step, starting from rest.
564      indic = 0
565      CALL dyn_spg( nit000, indic )       ! surface pressure gradient
566
567      ! the new velocity is ua*rdt
568
569      CALL lbc_lnk( ua, 'U', -1. )
570      CALL lbc_lnk( va, 'V', -1. )
571
572      ub(:,:,:) = ua(:,:,:) * rdt
573      vb(:,:,:) = va(:,:,:) * rdt
574      ua(:,:,:) = 0.e0
575      va(:,:,:) = 0.e0
576      un(:,:,:) = ub(:,:,:)
577      vn(:,:,:) = vb(:,:,:)
578       
579#if defined key_dynspg_rl
580      IF( lk_isl )   bsfb(:,:) = bsfn(:,:)          ! Put bsfb to zero
581#endif
582
583      ! Compute the divergence and curl
584
585      CALL div_cur( nit000 )            ! now horizontal divergence and curl
586
587      hdivb(:,:,:) = hdivn(:,:,:)       ! set the before to the now value
588      rotb (:,:,:) = rotn (:,:,:)       ! set the before to the now value
589      !
590   END SUBROUTINE istate_uvg
591
592   !!=====================================================================
593END MODULE istate
Note: See TracBrowser for help on using the repository browser.