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/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90 @ 5870

Last change on this file since 5870 was 5870, checked in by acc, 8 years ago

Branch 2015/dev_r5803_NOC_WAD. Merge in trunk changes from 5803 to 5869 in preparation for merge. Also tidied and reorganised some wetting and drying code. Renamed wadlmt.F90 to wetdry.F90. Wetting drying code changes restricted to domzgr.F90, domvvl.F90 nemogcm.F90 sshwzv.F90, dynspg_ts.F90, wetdry.F90 and dynhpg.F90. Code passes full SETTE tests with ln_wd=.false.. Still awaiting test case for checking with ln_wd=.false.
  • Property svn:keywords set to Id
File size: 26.3 KB
Line 
1MODULE istate
2   !!======================================================================
3   !!                     ***  MODULE  istate  ***
4   !! Ocean state   :  initial state setting
5   !!=====================================================================
6   !! History :  OPA  !  1989-12  (P. Andrich)  Original code
7   !!            5.0  !  1991-11  (G. Madec)  rewritting
8   !!            6.0  !  1996-01  (G. Madec)  terrain following coordinates
9   !!            8.0  !  2001-09  (M. Levy, M. Ben Jelloul)  istate_eel
10   !!            8.0  !  2001-09  (M. Levy, M. Ben Jelloul)  istate_uvg
11   !!   NEMO     1.0  !  2003-08  (G. Madec, C. Talandier)  F90: Free form, modules + EEL R5
12   !!             -   !  2004-05  (A. Koch-Larrouy)  istate_gyre
13   !!            2.0  !  2006-07  (S. Masson)  distributed restart using iom
14   !!            3.3  !  2010-10  (C. Ethe) merge TRC-TRA
15   !!            3.4  !  2011-04  (G. Madec) Merge of dtatem and dtasal & suppression of tb,tn/sb,sn
16   !!----------------------------------------------------------------------
17
18   !!----------------------------------------------------------------------
19   !!   istate_init   : initial state setting
20   !!   istate_tem    : analytical profile for initial Temperature
21   !!   istate_sal    : analytical profile for initial Salinity
22   !!   istate_eel    : initial state setting of EEL R5 configuration
23   !!   istate_gyre   : initial state setting of GYRE configuration
24   !!   istate_uvg    : initial velocity in geostropic balance
25   !!----------------------------------------------------------------------
26   USE oce             ! ocean dynamics and active tracers
27   USE dom_oce         ! ocean space and time domain
28   USE c1d             ! 1D vertical configuration
29   USE daymod          ! calendar
30   USE eosbn2          ! eq. of state, Brunt Vaisala frequency (eos     routine)
31   USE ldftra          ! lateral physics: ocean active tracers
32   USE zdf_oce         ! ocean vertical physics
33   USE phycst          ! physical constants
34   USE dtatsd          ! data temperature and salinity   (dta_tsd routine)
35   USE dtauvd          ! data: U & V current             (dta_uvd routine)
36   USE domvvl          ! varying vertical mesh
37   USE dynspg_oce      ! pressure gradient schemes
38   USE dynspg_flt      ! filtered free surface
39   USE sol_oce         ! ocean solver variables
40   !
41   USE in_out_manager  ! I/O manager
42   USE iom             ! I/O library
43   USE lib_mpp         ! MPP library
44   USE restart         ! restart
45   USE wrk_nemo        ! Memory allocation
46   USE timing          ! Timing
47
48   IMPLICIT NONE
49   PRIVATE
50
51   PUBLIC   istate_init   ! routine called by step.F90
52
53   !! * Substitutions
54#  include "domzgr_substitute.h90"
55#  include "vectopt_loop_substitute.h90"
56   !!----------------------------------------------------------------------
57   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
58   !! $Id$
59   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
60   !!----------------------------------------------------------------------
61CONTAINS
62
63   SUBROUTINE istate_init
64      !!----------------------------------------------------------------------
65      !!                   ***  ROUTINE istate_init  ***
66      !!
67      !! ** Purpose :   Initialization of the dynamics and tracer fields.
68      !!----------------------------------------------------------------------
69      INTEGER ::   ji, jj, jk   ! dummy loop indices
70      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zuvd    ! U & V data workspace
71      !!----------------------------------------------------------------------
72      !
73      IF( nn_timing == 1 )   CALL timing_start('istate_init')
74      !
75
76      IF(lwp) WRITE(numout,*)
77      IF(lwp) WRITE(numout,*) 'istate_ini : Initialization of the dynamics and tracers'
78      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
79
80                     CALL dta_tsd_init        ! Initialisation of T & S input data
81      IF( lk_c1d )   CALL dta_uvd_init        ! Initialization of U & V input data
82
83      rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk
84      rn2b (:,:,:  ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk
85      tsa  (:,:,:,:) = 0._wp                                   ! set one for all to 0 at level jpk
86      rab_b(:,:,:,:) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk
87
88      IF( ln_rstart ) THEN                    ! Restart from a file
89         !                                    ! -------------------
90         CALL rst_read                           ! Read the restart file
91         CALL day_init                           ! model calendar (using both namelist and restart infos)
92      ELSE
93         !                                    ! Start from rest
94         !                                    ! ---------------
95         numror = 0                              ! define numror = 0 -> no restart file to read
96         neuler = 0                              ! Set time-step indicator at nit000 (euler forward)
97         CALL day_init                           ! model calendar (using both namelist and restart infos)
98         !                                       ! Initialization of ocean to zero
99         !   before fields      !       now fields     
100         sshb (:,:)   = 0._wp   ;   sshn (:,:)   = 0._wp
101         ub   (:,:,:) = 0._wp   ;   un   (:,:,:) = 0._wp
102         vb   (:,:,:) = 0._wp   ;   vn   (:,:,:) = 0._wp 
103                                    hdivn(:,:,:) = 0._wp
104         !
105         IF( cp_cfg == 'eel' ) THEN
106            CALL istate_eel                      ! EEL   configuration : start from pre-defined U,V T-S fields
107         ELSEIF( cp_cfg == 'gyre' ) THEN         
108            CALL istate_gyre                     ! GYRE  configuration : start from pre-defined T-S fields
109         ELSE                                    ! Initial T-S, U-V fields read in files
110            IF ( ln_tsd_init ) THEN              ! read 3D T and S data at nit000
111               CALL dta_tsd( nit000, tsb ) 
112               tsn(:,:,:,:) = tsb(:,:,:,:)
113               !
114            ELSE                                 ! Initial T-S fields defined analytically
115               CALL istate_t_s
116            ENDIF
117            IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000
118               CALL wrk_alloc( jpi,jpj,jpk,2,   zuvd )
119               CALL dta_uvd( nit000, zuvd )
120               ub(:,:,:) = zuvd(:,:,:,1) ;  un(:,:,:) = ub(:,:,:)
121               vb(:,:,:) = zuvd(:,:,:,2) ;  vn(:,:,:) = vb(:,:,:)
122               CALL wrk_dealloc( jpi,jpj,jpk,2,   zuvd )
123            ENDIF
124         ENDIF
125         !   
126         ! - ML - sshn could be modified by istate_eel, so that initialization of fse3t_b is done here
127         IF( lk_vvl ) THEN
128            DO jk = 1, jpk
129               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
130            END DO
131         ENDIF
132         !
133      ENDIF
134      !
135      IF( lk_agrif ) THEN                  ! read free surface arrays in restart file
136         IF( ln_rstart ) THEN
137            IF( lk_dynspg_flt )  THEN      ! read or initialize the following fields
138               !                           ! gcx, gcxb for agrif_opa_init
139               IF( sol_oce_alloc()  > 0 )   CALL ctl_stop('agrif sol_oce_alloc: allocation of arrays failed')
140               CALL flt_rst( nit000, 'READ' )
141            ENDIF
142         ENDIF                             ! explicit case not coded yet with AGRIF
143      ENDIF
144      !
145      !
146      ! Initialize "now" and "before" barotropic velocities:
147      ! Do it whatever the free surface method, these arrays
148      ! being eventually used
149      !
150      !
151      un_b(:,:) = 0._wp   ;   vn_b(:,:) = 0._wp
152      ub_b(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp
153      !
154      DO jk = 1, jpkm1
155         DO jj = 1, jpj
156            DO ji = 1, jpi
157               un_b(ji,jj) = un_b(ji,jj) + fse3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
158               vn_b(ji,jj) = vn_b(ji,jj) + fse3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
159               !
160               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
161               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
162            END DO
163         END DO
164      END DO
165      !
166      un_b(:,:) = un_b(:,:) * hur  (:,:)
167      vn_b(:,:) = vn_b(:,:) * hvr  (:,:)
168      !
169      ub_b(:,:) = ub_b(:,:) * hur_b(:,:)
170      vb_b(:,:) = vb_b(:,:) * hvr_b(:,:)
171      !
172      !
173      IF( nn_timing == 1 )   CALL timing_stop('istate_init')
174      !
175   END SUBROUTINE istate_init
176
177
178   SUBROUTINE istate_t_s
179      !!---------------------------------------------------------------------
180      !!                  ***  ROUTINE istate_t_s  ***
181      !!   
182      !! ** Purpose :   Intialization of the temperature field with an
183      !!      analytical profile or a file (i.e. in EEL configuration)
184      !!
185      !! ** Method  : - temperature: use Philander analytic profile
186      !!              - salinity   : use to a constant value 35.5
187      !!
188      !! References :  Philander ???
189      !!----------------------------------------------------------------------
190      INTEGER  ::   ji, jj, jk
191      REAL(wp) ::   zsal = 35.50_wp
192      !!----------------------------------------------------------------------
193      !
194      IF(lwp) WRITE(numout,*)
195      IF(lwp) WRITE(numout,*) 'istate_t_s : Philander s initial temperature profile'
196      IF(lwp) WRITE(numout,*) '~~~~~~~~~~   and constant salinity (',zsal,' psu)'
197      !
198      DO jk = 1, jpk
199         tsn(:,:,jk,jp_tem) = (  ( ( 7.5 - 0. * ABS( gphit(:,:) )/30. ) * ( 1.-TANH((fsdept(:,:,jk)-80.)/30.) )   &
200            &                + 10. * ( 5000. - fsdept(:,:,jk) ) /5000.)  ) * tmask(:,:,jk)
201         tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem)
202      END DO
203      tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)
204      tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
205      !
206   END SUBROUTINE istate_t_s
207
208
209   SUBROUTINE istate_eel
210      !!----------------------------------------------------------------------
211      !!                   ***  ROUTINE istate_eel  ***
212      !!
213      !! ** Purpose :   Initialization of the dynamics and tracers for EEL R5
214      !!      configuration (channel with or without a topographic bump)
215      !!
216      !! ** Method  : - set temprature field
217      !!              - set salinity field
218      !!              - set velocity field including horizontal divergence
219      !!                and relative vorticity fields
220      !!----------------------------------------------------------------------
221      USE divhor     ! hor. divergence      (div_hor routine)
222      USE iom
223      !
224      INTEGER  ::   inum              ! temporary logical unit
225      INTEGER  ::   ji, jj, jk        ! dummy loop indices
226      INTEGER  ::   ijloc
227      REAL(wp) ::   zh1, zh2, zslope, zcst, zfcor   ! temporary scalars
228      REAL(wp) ::   zt1  = 15._wp                   ! surface temperature value (EEL R5)
229      REAL(wp) ::   zt2  =  5._wp                   ! bottom  temperature value (EEL R5)
230      REAL(wp) ::   zsal = 35.0_wp                  ! constant salinity (EEL R2, R5 and R6)
231      REAL(wp) ::   zueel = 0.1_wp                  ! constant uniform zonal velocity (EEL R5)
232      REAL(wp), DIMENSION(jpiglo,jpjglo) ::   zssh  ! initial ssh over the global domain
233      !!----------------------------------------------------------------------
234      !
235      SELECT CASE ( jp_cfg ) 
236         !                                              ! ====================
237         CASE ( 5 )                                     ! EEL R5 configuration
238            !                                           ! ====================
239            !
240            ! set temperature field with a linear profile
241            ! -------------------------------------------
242            IF(lwp) WRITE(numout,*)
243            IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: linear temperature profile'
244            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
245            !
246            zh1 = gdept_1d(  1  )
247            zh2 = gdept_1d(jpkm1)
248            !
249            zslope = ( zt1 - zt2 ) / ( zh1 - zh2 )
250            zcst   = ( zt1 * ( zh1 - zh2) - ( zt1 - zt2 ) * zh1 ) / ( zh1 - zh2 )
251            !
252            DO jk = 1, jpk
253               tsn(:,:,jk,jp_tem) = ( zt2 + zt1 * exp( - fsdept(:,:,jk) / 1000 ) ) * tmask(:,:,jk)
254               tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem)
255            END DO
256            !
257            IF(lwp) CALL prizre( tsn(:,:,:,jp_tem), jpi   , jpj   , jpk   , jpj/2 ,   &
258               &                             1     , jpi   , 5     , 1     , jpk   ,   &
259               &                             1     , 1.    , numout                  )
260            !
261            ! set salinity field to a constant value
262            ! --------------------------------------
263            IF(lwp) WRITE(numout,*)
264            IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: constant salinity field, S = ', zsal
265            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
266            !
267            tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)
268            tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
269            !
270            ! set the dynamics: U,V, hdiv (and ssh if necessary)
271            ! ----------------
272            ! Start EEL5 configuration with barotropic geostrophic velocities
273            ! according the sshb and sshn SSH imposed.
274            ! we assume a uniform grid (hence the use of e1t(1,1) for delta_y)
275            ! we use the Coriolis frequency at mid-channel.   
276            ub(:,:,:) = zueel * umask(:,:,:)
277            un(:,:,:) = ub(:,:,:)
278            ijloc = mj0(INT(jpjglo-1)/2)
279            zfcor = ff(1,ijloc)
280            !
281            DO jj = 1, jpjglo
282               zssh(:,jj) = - (FLOAT(jj)- FLOAT(jpjglo-1)/2.)*zueel*e1t(1,1)*zfcor/grav 
283            END DO
284            !
285            IF(lwp) THEN
286               WRITE(numout,*) ' Uniform zonal velocity for EEL R5:',zueel
287               WRITE(numout,*) ' Geostrophic SSH profile as a function of y:'
288               WRITE(numout,'(12(1x,f6.2))') zssh(1,:)
289            ENDIF
290            !
291            DO jj = 1, nlcj
292               DO ji = 1, nlci
293                  sshb(ji,jj) = zssh( mig(ji) , mjg(jj) ) * tmask(ji,jj,1)
294               END DO
295            END DO
296            sshb(nlci+1:jpi,      :   ) = 0.e0      ! set to zero extra mpp columns
297            sshb(      :   ,nlcj+1:jpj) = 0.e0      ! set to zero extra mpp rows
298            !
299            sshn(:,:) = sshb(:,:)                   ! set now ssh to the before value
300            !
301            IF( nn_rstssh /= 0 ) THEN 
302               nn_rstssh = 0                        ! hand-made initilization of ssh
303               CALL ctl_warn( 'istate_eel: force nn_rstssh = 0' )
304            ENDIF
305            !
306!!gm  Check  here call to div_hor should not be necessary
307!!gm         div_hor call runoffs  not sure they are defined at that level
308            CALL div_hor( nit000 )                  ! horizontal divergence and relative vorticity (curl)
309            ! N.B. the vertical velocity will be computed from the horizontal divergence field
310            ! in istate by a call to wzv routine
311
312
313            !                                     ! ==========================
314         CASE ( 2 , 6 )                           ! EEL R2 or R6 configuration
315            !                                     ! ==========================
316            !
317            ! set temperature field with a NetCDF file
318            ! ----------------------------------------
319            IF(lwp) WRITE(numout,*)
320            IF(lwp) WRITE(numout,*) 'istate_eel : EEL R2 or R6: read initial temperature in a NetCDF file'
321            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
322            !
323            CALL iom_open ( 'eel.initemp', inum )
324            CALL iom_get ( inum, jpdom_data, 'initemp', tsb(:,:,:,jp_tem) ) ! read before temprature (tb)
325            CALL iom_close( inum )
326            !
327            tsn(:,:,:,jp_tem) = tsb(:,:,:,jp_tem)                            ! set nox temperature to tb
328            !
329            IF(lwp) CALL prizre( tsn(:,:,:,jp_tem), jpi   , jpj   , jpk   , jpj/2 ,   &
330               &                            1     , jpi   , 5     , 1     , jpk   ,   &
331               &                            1     , 1.    , numout                  )
332            !
333            ! set salinity field to a constant value
334            ! --------------------------------------
335            IF(lwp) WRITE(numout,*)
336            IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: constant salinity field, S = ', zsal
337            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
338            !
339            tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)
340            tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
341            !
342            !                                    ! ===========================
343         CASE DEFAULT                            ! NONE existing configuration
344            !                                    ! ===========================
345            WRITE(ctmp1,*) 'EEL with a ', jp_cfg,' km resolution is not coded'
346            CALL ctl_stop( ctmp1 )
347            !
348      END SELECT
349      !
350   END SUBROUTINE istate_eel
351
352
353   SUBROUTINE istate_gyre
354      !!----------------------------------------------------------------------
355      !!                   ***  ROUTINE istate_gyre  ***
356      !!
357      !! ** Purpose :   Initialization of the dynamics and tracers for GYRE
358      !!      configuration (double gyre with rotated domain)
359      !!
360      !! ** Method  : - set temprature field
361      !!              - set salinity   field
362      !!----------------------------------------------------------------------
363      INTEGER :: ji, jj, jk  ! dummy loop indices
364      INTEGER            ::   inum          ! temporary logical unit
365      INTEGER, PARAMETER ::   ntsinit = 0   ! (0/1) (analytical/input data files) T&S initialization
366      !!----------------------------------------------------------------------
367      !
368      SELECT CASE ( ntsinit)
369      !
370      CASE ( 0 )                  ! analytical T/S profil deduced from LEVITUS
371         IF(lwp) WRITE(numout,*)
372         IF(lwp) WRITE(numout,*) 'istate_gyre : initial analytical T and S profil deduced from LEVITUS '
373         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
374         !
375         DO jk = 1, jpk
376            DO jj = 1, jpj
377               DO ji = 1, jpi
378                  tsn(ji,jj,jk,jp_tem) = (  16. - 12. * TANH( (fsdept(ji,jj,jk) - 400) / 700 )         )   &
379                       &           * (-TANH( (500-fsdept(ji,jj,jk)) / 150 ) + 1) / 2               &
380                       &       + (      15. * ( 1. - TANH( (fsdept(ji,jj,jk)-50.) / 1500.) )       &
381                       &                - 1.4 * TANH((fsdept(ji,jj,jk)-100.) / 100.)               &   
382                       &                + 7.  * (1500. - fsdept(ji,jj,jk)) / 1500.             )   & 
383                       &           * (-TANH( (fsdept(ji,jj,jk) - 500) / 150) + 1) / 2
384                  tsn(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
385                  tsb(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem)
386
387                  tsn(ji,jj,jk,jp_sal) =  (  36.25 - 1.13 * TANH( (fsdept(ji,jj,jk) - 305) / 460 )  )  &
388                     &              * (-TANH((500 - fsdept(ji,jj,jk)) / 150) + 1) / 2          &
389                     &          + (  35.55 + 1.25 * (5000. - fsdept(ji,jj,jk)) / 5000.         &
390                     &                - 1.62 * TANH( (fsdept(ji,jj,jk) - 60.  ) / 650. )       &
391                     &                + 0.2  * TANH( (fsdept(ji,jj,jk) - 35.  ) / 100. )       &
392                     &                + 0.2  * TANH( (fsdept(ji,jj,jk) - 1000.) / 5000.)    )  &
393                     &              * (-TANH((fsdept(ji,jj,jk) - 500) / 150) + 1) / 2 
394                  tsn(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
395                  tsb(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal)
396               END DO
397            END DO
398         END DO
399         !
400      CASE ( 1 )                  ! T/S data fields read in dta_tem.nc/data_sal.nc files
401         IF(lwp) WRITE(numout,*)
402         IF(lwp) WRITE(numout,*) 'istate_gyre : initial T and S read from dta_tem.nc/data_sal.nc files'
403         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
404         IF(lwp) WRITE(numout,*) '              NetCDF FORMAT'
405
406         ! Read temperature field
407         ! ----------------------
408         CALL iom_open ( 'data_tem', inum )
409         CALL iom_get ( inum, jpdom_data, 'votemper', tsn(:,:,:,jp_tem) ) 
410         CALL iom_close( inum )
411
412         tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:) 
413         tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
414
415         ! Read salinity field
416         ! -------------------
417         CALL iom_open ( 'data_sal', inum )
418         CALL iom_get ( inum, jpdom_data, 'vosaline', tsn(:,:,:,jp_sal) ) 
419         CALL iom_close( inum )
420
421         tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:) 
422         tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
423         !
424      END SELECT
425      !
426      IF(lwp) THEN
427         WRITE(numout,*)
428         WRITE(numout,*) '              Initial temperature and salinity profiles:'
429         WRITE(numout, "(9x,' level   gdept_1d   temperature   salinity   ')" )
430         WRITE(numout, "(10x, i4, 3f10.2)" ) ( jk, gdept_1d(jk), tsn(2,2,jk,jp_tem), tsn(2,2,jk,jp_sal), jk = 1, jpk )
431      ENDIF
432      !
433   END SUBROUTINE istate_gyre
434
435
436   SUBROUTINE istate_uvg
437      !!----------------------------------------------------------------------
438      !!                  ***  ROUTINE istate_uvg  ***
439      !!
440      !! ** Purpose :   Compute the geostrophic velocities from (tn,sn) fields
441      !!
442      !! ** Method  :   Using the hydrostatic hypothesis the now hydrostatic
443      !!      pressure is computed by integrating the in-situ density from the
444      !!      surface to the bottom.
445      !!                 p=integral [ rau*g dz ]
446      !!----------------------------------------------------------------------
447      USE dynspg          ! surface pressure gradient             (dyn_spg routine)
448      USE divhor          ! hor. divergence                       (div_hor routine)
449      USE lbclnk          ! ocean lateral boundary condition (or mpp link)
450      !
451      INTEGER ::   ji, jj, jk        ! dummy loop indices
452      INTEGER ::   indic             ! ???
453      REAL(wp) ::   zmsv, zphv, zmsu, zphu, zalfg     ! temporary scalars
454      REAL(wp), POINTER, DIMENSION(:,:,:) :: zprn
455      !!----------------------------------------------------------------------
456      !
457      CALL wrk_alloc( jpi,jpj,jpk,   zprn)
458      !
459      IF(lwp) WRITE(numout,*) 
460      IF(lwp) WRITE(numout,*) 'istate_uvg : Start from Geostrophy'
461      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
462
463      ! Compute the now hydrostatic pressure
464      ! ------------------------------------
465
466      zalfg = 0.5 * grav * rau0
467     
468      zprn(:,:,1) = zalfg * fse3w(:,:,1) * ( 1 + rhd(:,:,1) )       ! Surface value
469
470      DO jk = 2, jpkm1                                              ! Vertical integration from the surface
471         zprn(:,:,jk) = zprn(:,:,jk-1)   &
472            &         + zalfg * fse3w(:,:,jk) * ( 2. + rhd(:,:,jk) + rhd(:,:,jk-1) )
473      END DO 
474
475      ! Compute geostrophic balance
476      ! ---------------------------
477      DO jk = 1, jpkm1
478         DO jj = 2, jpjm1
479            DO ji = fs_2, fs_jpim1   ! vertor opt.
480               zmsv = 1. / MAX(  umask(ji-1,jj+1,jk) + umask(ji  ,jj+1,jk)   &
481                               + umask(ji-1,jj  ,jk) + umask(ji  ,jj  ,jk) , 1.  )
482               zphv = ( zprn(ji  ,jj+1,jk) - zprn(ji-1,jj+1,jk) ) * umask(ji-1,jj+1,jk) / e1u(ji-1,jj+1)   &
483                    + ( zprn(ji+1,jj+1,jk) - zprn(ji  ,jj+1,jk) ) * umask(ji  ,jj+1,jk) / e1u(ji  ,jj+1)   &
484                    + ( zprn(ji  ,jj  ,jk) - zprn(ji-1,jj  ,jk) ) * umask(ji-1,jj  ,jk) / e1u(ji-1,jj  )   &
485                    + ( zprn(ji+1,jj  ,jk) - zprn(ji  ,jj  ,jk) ) * umask(ji  ,jj  ,jk) / e1u(ji  ,jj  )
486               zphv = 1. / rau0 * zphv * zmsv * vmask(ji,jj,jk)
487
488               zmsu = 1. / MAX(  vmask(ji+1,jj  ,jk) + vmask(ji  ,jj  ,jk)   &
489                               + vmask(ji+1,jj-1,jk) + vmask(ji  ,jj-1,jk) , 1.  )
490               zphu = ( zprn(ji+1,jj+1,jk) - zprn(ji+1,jj  ,jk) ) * vmask(ji+1,jj  ,jk) / e2v(ji+1,jj  )   &
491                    + ( zprn(ji  ,jj+1,jk) - zprn(ji  ,jj  ,jk) ) * vmask(ji  ,jj  ,jk) / e2v(ji  ,jj  )   &
492                    + ( zprn(ji+1,jj  ,jk) - zprn(ji+1,jj-1,jk) ) * vmask(ji+1,jj-1,jk) / e2v(ji+1,jj-1)   &
493                    + ( zprn(ji  ,jj  ,jk) - zprn(ji  ,jj-1,jk) ) * vmask(ji  ,jj-1,jk) / e2v(ji  ,jj-1)
494               zphu = 1. / rau0 * zphu * zmsu * umask(ji,jj,jk)
495
496               ! Compute the geostrophic velocities
497               un(ji,jj,jk) = -2. * zphu / ( ff(ji,jj) + ff(ji  ,jj-1) )
498               vn(ji,jj,jk) =  2. * zphv / ( ff(ji,jj) + ff(ji-1,jj  ) )
499            END DO
500         END DO
501      END DO
502
503      IF(lwp) WRITE(numout,*) '         we force to zero bottom velocity'
504
505      ! Susbtract the bottom velocity (level jpk-1 for flat bottom case)
506      ! to have a zero bottom velocity
507
508      DO jk = 1, jpkm1
509         un(:,:,jk) = ( un(:,:,jk) - un(:,:,jpkm1) ) * umask(:,:,jk)
510         vn(:,:,jk) = ( vn(:,:,jk) - vn(:,:,jpkm1) ) * vmask(:,:,jk)
511      END DO
512
513      CALL lbc_lnk( un, 'U', -1. )
514      CALL lbc_lnk( vn, 'V', -1. )
515     
516      ub(:,:,:) = un(:,:,:)
517      vb(:,:,:) = vn(:,:,:)
518     
519      ! WARNING !!!!!
520      ! after initializing u and v, we need to calculate the initial streamfunction bsf.
521      ! Otherwise, only the trend will be computed and the model will blow up (inconsistency).
522      ! to do that, we call dyn_spg with a special trick:
523      ! we fill ua and va with the velocities divided by dt, and the streamfunction will be brought to the
524      ! right value assuming the velocities have been set up in one time step.
525      ! we then set bsfd to zero (first guess for next step is d(psi)/dt = 0.)
526      !  sets up s false trend to calculate the barotropic streamfunction.
527
528      ua(:,:,:) = ub(:,:,:) / rdt
529      va(:,:,:) = vb(:,:,:) / rdt
530
531      ! calls dyn_spg. we assume euler time step, starting from rest.
532      indic = 0
533      CALL dyn_spg( nit000, indic )       ! surface pressure gradient
534      !
535      ! the new velocity is ua*rdt
536      !
537      CALL lbc_lnk( ua, 'U', -1. )
538      CALL lbc_lnk( va, 'V', -1. )
539
540      ub(:,:,:) = ua(:,:,:) * rdt
541      vb(:,:,:) = va(:,:,:) * rdt
542      ua(:,:,:) = 0.e0
543      va(:,:,:) = 0.e0
544      un(:,:,:) = ub(:,:,:)
545      vn(:,:,:) = vb(:,:,:)
546      !
547!!gm  Check  here call to div_hor should not be necessary
548!!gm         div_hor call runoffs  not sure they are defined at that level
549      CALL div_hor( nit000 )            ! now horizontal divergence
550      !
551      CALL wrk_dealloc( jpi,jpj,jpk,   zprn)
552      !
553   END SUBROUTINE istate_uvg
554
555   !!=====================================================================
556END MODULE istate
Note: See TracBrowser for help on using the repository browser.