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/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90 @ 3211

Last change on this file since 3211 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

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