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 trunk/NEMO/OPA_SRC – NEMO

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

Last change on this file since 1200 was 1200, checked in by rblod, 16 years ago

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