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.
istate.F90 in branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC – NEMO

source: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/istate.F90 @ 2066

Last change on this file since 2066 was 2066, checked in by rblod, 14 years ago

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