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/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90 @ 11101

Last change on this file since 11101 was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.
File size: 27.6 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_oce      ! ocean active tracers: lateral physics
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 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      ! filtered free surface
41   USE sol_oce         ! ocean solver variables
42   !
43   USE in_out_manager  ! I/O manager
44   USE iom             ! I/O library
45   USE lib_mpp         ! MPP library
46   USE restart         ! restart
47   USE wrk_nemo        ! Memory allocation
48   USE timing          ! Timing
49
50   IMPLICIT NONE
51   PRIVATE
52
53   PUBLIC   istate_init   ! routine called by step.F90
54
55   !! * Substitutions
56#  include "domzgr_substitute.h90"
57#  include "vectopt_loop_substitute.h90"
58   !!----------------------------------------------------------------------
59   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
60   !! $Id$
61   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
62   !!----------------------------------------------------------------------
63CONTAINS
64
65   SUBROUTINE istate_init
66      !!----------------------------------------------------------------------
67      !!                   ***  ROUTINE istate_init  ***
68      !!
69      !! ** Purpose :   Initialization of the dynamics and tracer fields.
70      !!----------------------------------------------------------------------
71      INTEGER ::   ji, jj, jk   ! dummy loop indices
72      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zuvd    ! U & V data workspace
73      !!----------------------------------------------------------------------
74      !
75      IF( nn_timing == 1 )   CALL timing_start('istate_init')
76      !
77
78      IF(lwp) WRITE(numout,*) ' '
79      IF(lwp) WRITE(numout,*) 'istate_ini : Initialization of the dynamics and tracers'
80      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
81      IF(lwp .AND. lflush) CALL flush(numout)
82
83      CALL dta_tsd_init                       ! Initialisation of T & S input data
84      IF( lk_c1d ) CALL dta_uvd_init          ! Initialization of U & V input data
85
86      rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk
87      rn2b (:,:,:  ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk
88      tsa  (:,:,:,:) = 0._wp                                   ! set one for all to 0 at level jpk
89      rab_b(:,:,:,:) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk
90
91      IF( ln_rstart ) THEN                    ! Restart from a file
92         !                                    ! -------------------
93         CALL rst_read                           ! Read the restart file
94         CALL day_init                           ! model calendar (using both namelist and restart infos)
95      ELSE
96         !                                    ! Start from rest
97         !                                    ! ---------------
98         numror = 0                              ! define numror = 0 -> no restart file to read
99         neuler = 0                              ! Set time-step indicator at nit000 (euler forward)
100         CALL day_init                           ! model calendar (using both namelist and restart infos)
101         !                                       ! Initialization of ocean to zero
102         !   before fields      !       now fields     
103         sshb (:,:)   = 0._wp   ;   sshn (:,:)   = 0._wp
104         ub   (:,:,:) = 0._wp   ;   un   (:,:,:) = 0._wp
105         vb   (:,:,:) = 0._wp   ;   vn   (:,:,:) = 0._wp 
106         rotb (:,:,:) = 0._wp   ;   rotn (:,:,:) = 0._wp
107         hdivb(:,:,:) = 0._wp   ;   hdivn(:,:,:) = 0._wp
108         !
109         IF( cp_cfg == 'eel' ) THEN
110            CALL istate_eel                      ! EEL   configuration : start from pre-defined U,V T-S fields
111         ELSEIF( cp_cfg == 'gyre' ) THEN         
112            CALL istate_gyre                     ! GYRE  configuration : start from pre-defined T-S fields
113         ELSE                                    ! Initial T-S, U-V fields read in files
114            IF ( ln_tsd_init ) THEN              ! read 3D T and S data at nit000
115               CALL dta_tsd( nit000, tsb ) 
116               tsn(:,:,:,:) = tsb(:,:,:,:)
117               !
118            ELSE                                 ! Initial T-S fields defined analytically
119               CALL istate_t_s
120            ENDIF
121            IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000
122               CALL wrk_alloc( jpi, jpj, jpk, 2, zuvd )
123               CALL dta_uvd( nit000, zuvd )
124               ub(:,:,:) = zuvd(:,:,:,1) ;  un(:,:,:) = ub(:,:,:)
125               vb(:,:,:) = zuvd(:,:,:,2) ;  vn(:,:,:) = vb(:,:,:)
126               CALL wrk_dealloc( jpi, jpj, jpk, 2, zuvd )
127            ENDIF
128         ENDIF
129         !
130         CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) )        ! before potential and in situ densities
131#if ! defined key_c1d
132         IF( ln_zps .AND. .NOT. ln_isfcav)                                 &
133            &            CALL zps_hde    ( nit000, jpts, tsb, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
134            &                                            rhd, gru , grv    )  ! of t, s, rd at the last ocean level
135         IF( ln_zps .AND.       ln_isfcav)                                 &
136            &            CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
137            &                                            rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
138            &                                     gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
139#endif
140         !   
141         ! - ML - sshn could be modified by istate_eel, so that initialization of fse3t_b is done here
142         IF( lk_vvl ) THEN
143            DO jk = 1, jpk
144               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
145            ENDDO
146         ENDIF
147         !
148      ENDIF
149      !
150      IF( lk_agrif ) THEN                  ! read free surface arrays in restart file
151         IF( ln_rstart ) THEN
152            IF( lk_dynspg_flt )  THEN      ! read or initialize the following fields
153               !                           ! gcx, gcxb for agrif_opa_init
154               IF( sol_oce_alloc()  > 0 )   CALL ctl_stop('agrif sol_oce_alloc: allocation of arrays failed')
155               CALL flt_rst( nit000, 'READ' )
156            ENDIF
157         ENDIF                             ! explicit case not coded yet with AGRIF
158      ENDIF
159      !
160      !
161      ! Initialize "now" and "before" barotropic velocities:
162      ! Do it whatever the free surface method, these arrays
163      ! being eventually used
164      !
165      !
166      un_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp
167      ub_b(:,:) = 0._wp ; vb_b(:,:) = 0._wp
168      !
169      DO jk = 1, jpkm1
170         DO jj = 1, jpj
171            DO ji = 1, jpi
172               un_b(ji,jj) = un_b(ji,jj) + fse3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
173               vn_b(ji,jj) = vn_b(ji,jj) + fse3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
174               !
175               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
176               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
177            END DO
178         END DO
179      END DO
180      !
181      un_b(:,:) = un_b(:,:) * hur  (:,:)
182      vn_b(:,:) = vn_b(:,:) * hvr  (:,:)
183      !
184      ub_b(:,:) = ub_b(:,:) * hur_b(:,:)
185      vb_b(:,:) = vb_b(:,:) * hvr_b(:,:)
186      !
187      !
188      IF( nn_timing == 1 )   CALL timing_stop('istate_init')
189      !
190   END SUBROUTINE istate_init
191
192
193   SUBROUTINE istate_t_s
194      !!---------------------------------------------------------------------
195      !!                  ***  ROUTINE istate_t_s  ***
196      !!   
197      !! ** Purpose :   Intialization of the temperature field with an
198      !!      analytical profile or a file (i.e. in EEL configuration)
199      !!
200      !! ** Method  : - temperature: use Philander analytic profile
201      !!              - salinity   : use to a constant value 35.5
202      !!
203      !! References :  Philander ???
204      !!----------------------------------------------------------------------
205      INTEGER  :: ji, jj, jk
206      REAL(wp) ::   zsal = 35.50
207      !!----------------------------------------------------------------------
208      !
209      IF(lwp) WRITE(numout,*)
210      IF(lwp) WRITE(numout,*) 'istate_t_s : Philander s initial temperature profile'
211      IF(lwp) WRITE(numout,*) '~~~~~~~~~~   and constant salinity (',zsal,' psu)'
212      IF(lwp .AND. lflush) CALL flush(numout)
213      !
214      DO jk = 1, jpk
215         tsn(:,:,jk,jp_tem) = (  ( ( 7.5 - 0. * ABS( gphit(:,:) )/30. ) * ( 1.-TANH((fsdept(:,:,jk)-80.)/30.) )   &
216            &                + 10. * ( 5000. - fsdept(:,:,jk) ) /5000.)  ) * tmask(:,:,jk)
217         tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem)
218      END DO
219      tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)
220      tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
221      !
222   END SUBROUTINE istate_t_s
223
224
225   SUBROUTINE istate_eel
226      !!----------------------------------------------------------------------
227      !!                   ***  ROUTINE istate_eel  ***
228      !!
229      !! ** Purpose :   Initialization of the dynamics and tracers for EEL R5
230      !!      configuration (channel with or without a topographic bump)
231      !!
232      !! ** Method  : - set temprature field
233      !!              - set salinity field
234      !!              - set velocity field including horizontal divergence
235      !!                and relative vorticity fields
236      !!----------------------------------------------------------------------
237      USE divcur     ! hor. divergence & rel. vorticity      (div_cur routine)
238      USE iom
239      !
240      INTEGER  ::   inum              ! temporary logical unit
241      INTEGER  ::   ji, jj, jk        ! dummy loop indices
242      INTEGER  ::   ijloc
243      REAL(wp) ::   zh1, zh2, zslope, zcst, zfcor   ! temporary scalars
244      REAL(wp) ::   zt1  = 15._wp                   ! surface temperature value (EEL R5)
245      REAL(wp) ::   zt2  =  5._wp                   ! bottom  temperature value (EEL R5)
246      REAL(wp) ::   zsal = 35.0_wp                  ! constant salinity (EEL R2, R5 and R6)
247      REAL(wp) ::   zueel = 0.1_wp                  ! constant uniform zonal velocity (EEL R5)
248      REAL(wp), DIMENSION(jpiglo,jpjglo) ::   zssh  ! initial ssh over the global domain
249      !!----------------------------------------------------------------------
250      !
251      SELECT CASE ( jp_cfg ) 
252         !                                              ! ====================
253         CASE ( 5 )                                     ! EEL R5 configuration
254            !                                           ! ====================
255            !
256            ! set temperature field with a linear profile
257            ! -------------------------------------------
258            IF(lwp) WRITE(numout,*)
259            IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: linear temperature profile'
260            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
261            !
262            zh1 = gdept_1d(  1  )
263            zh2 = gdept_1d(jpkm1)
264            !
265            zslope = ( zt1 - zt2 ) / ( zh1 - zh2 )
266            zcst   = ( zt1 * ( zh1 - zh2) - ( zt1 - zt2 ) * zh1 ) / ( zh1 - zh2 )
267            !
268            DO jk = 1, jpk
269               tsn(:,:,jk,jp_tem) = ( zt2 + zt1 * exp( - fsdept(:,:,jk) / 1000 ) ) * tmask(:,:,jk)
270               tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem)
271            END DO
272            !
273            IF(lwp) CALL prizre( tsn(:,:,:,jp_tem), jpi   , jpj   , jpk   , jpj/2 ,   &
274               &                             1     , jpi   , 5     , 1     , jpk   ,   &
275               &                             1     , 1.    , numout                  )
276            !
277            ! set salinity field to a constant value
278            ! --------------------------------------
279            IF(lwp) WRITE(numout,*)
280            IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: constant salinity field, S = ', zsal
281            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
282            !
283            tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)
284            tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
285            !
286            ! set the dynamics: U,V, hdiv, rot (and ssh if necessary)
287            ! ----------------
288            ! Start EEL5 configuration with barotropic geostrophic velocities
289            ! according the sshb and sshn SSH imposed.
290            ! we assume a uniform grid (hence the use of e1t(1,1) for delta_y)
291            ! we use the Coriolis frequency at mid-channel.   
292            ub(:,:,:) = zueel * umask(:,:,:)
293            un(:,:,:) = ub(:,:,:)
294            ijloc = mj0(INT(jpjglo-1)/2)
295            zfcor = ff(1,ijloc)
296            !
297            DO jj = 1, jpjglo
298               zssh(:,jj) = - (FLOAT(jj)- FLOAT(jpjglo-1)/2.)*zueel*e1t(1,1)*zfcor/grav 
299            END DO
300            !
301            IF(lwp) THEN
302               WRITE(numout,*) ' Uniform zonal velocity for EEL R5:',zueel
303               WRITE(numout,*) ' Geostrophic SSH profile as a function of y:'
304               WRITE(numout,'(12(1x,f6.2))') zssh(1,:)
305            ENDIF
306            !
307            DO jj = 1, nlcj
308               DO ji = 1, nlci
309                  sshb(ji,jj) = zssh( mig(ji) , mjg(jj) ) * tmask(ji,jj,1)
310               END DO
311            END DO
312            sshb(nlci+1:jpi,      :   ) = 0.e0      ! set to zero extra mpp columns
313            sshb(      :   ,nlcj+1:jpj) = 0.e0      ! set to zero extra mpp rows
314            !
315            sshn(:,:) = sshb(:,:)                   ! set now ssh to the before value
316            !
317            IF( nn_rstssh /= 0 ) THEN 
318               nn_rstssh = 0                        ! hand-made initilization of ssh
319               CALL ctl_warn( 'istate_eel: force nn_rstssh = 0' )
320            ENDIF
321            !
322            CALL div_cur( nit000 )                  ! horizontal divergence and relative vorticity (curl)
323            ! N.B. the vertical velocity will be computed from the horizontal divergence field
324            ! in istate by a call to wzv routine
325
326
327            !                                     ! ==========================
328         CASE ( 2 , 6 )                           ! EEL R2 or R6 configuration
329            !                                     ! ==========================
330            !
331            ! set temperature field with a NetCDF file
332            ! ----------------------------------------
333            IF(lwp) WRITE(numout,*)
334            IF(lwp) WRITE(numout,*) 'istate_eel : EEL R2 or R6: read initial temperature in a NetCDF file'
335            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
336            !
337            CALL iom_open ( 'eel.initemp', inum )
338            CALL iom_get ( inum, jpdom_data, 'initemp', tsb(:,:,:,jp_tem) ) ! read before temprature (tb)
339            CALL iom_close( inum )
340            !
341            tsn(:,:,:,jp_tem) = tsb(:,:,:,jp_tem)                            ! set nox temperature to tb
342            !
343            IF(lwp) CALL prizre( tsn(:,:,:,jp_tem), jpi   , jpj   , jpk   , jpj/2 ,   &
344               &                            1     , jpi   , 5     , 1     , jpk   ,   &
345               &                            1     , 1.    , numout                  )
346            !
347            ! set salinity field to a constant value
348            ! --------------------------------------
349            IF(lwp) WRITE(numout,*)
350            IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: constant salinity field, S = ', zsal
351            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
352            !
353            tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)
354            tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
355            !
356            !                                    ! ===========================
357         CASE DEFAULT                            ! NONE existing configuration
358            !                                    ! ===========================
359            WRITE(ctmp1,*) 'EEL with a ', jp_cfg,' km resolution is not coded'
360            CALL ctl_stop( ctmp1 )
361            !
362      END SELECT
363      !
364      IF(lwp .AND. lflush) CALL flush(numout)
365      !
366   END SUBROUTINE istate_eel
367
368
369   SUBROUTINE istate_gyre
370      !!----------------------------------------------------------------------
371      !!                   ***  ROUTINE istate_gyre  ***
372      !!
373      !! ** Purpose :   Initialization of the dynamics and tracers for GYRE
374      !!      configuration (double gyre with rotated domain)
375      !!
376      !! ** Method  : - set temprature field
377      !!              - set salinity field
378      !!----------------------------------------------------------------------
379      INTEGER :: ji, jj, jk  ! dummy loop indices
380      INTEGER            ::   inum          ! temporary logical unit
381      INTEGER, PARAMETER ::   ntsinit = 0   ! (0/1) (analytical/input data files) T&S initialization
382      !!----------------------------------------------------------------------
383      !
384      SELECT CASE ( ntsinit)
385      !
386      CASE ( 0 )                  ! analytical T/S profil deduced from LEVITUS
387         IF(lwp) WRITE(numout,*)
388         IF(lwp) WRITE(numout,*) 'istate_gyre : initial analytical T and S profil deduced from LEVITUS '
389         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
390         !
391         DO jk = 1, jpk
392            DO jj = 1, jpj
393               DO ji = 1, jpi
394                  tsn(ji,jj,jk,jp_tem) = (  16. - 12. * TANH( (fsdept(ji,jj,jk) - 400) / 700 )         )   &
395                       &           * (-TANH( (500-fsdept(ji,jj,jk)) / 150 ) + 1) / 2               &
396                       &       + (      15. * ( 1. - TANH( (fsdept(ji,jj,jk)-50.) / 1500.) )       &
397                       &                - 1.4 * TANH((fsdept(ji,jj,jk)-100.) / 100.)               &   
398                       &                + 7.  * (1500. - fsdept(ji,jj,jk)) / 1500.             )   & 
399                       &           * (-TANH( (fsdept(ji,jj,jk) - 500) / 150) + 1) / 2
400                  tsn(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
401                  tsb(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem)
402
403                  tsn(ji,jj,jk,jp_sal) =  (  36.25 - 1.13 * TANH( (fsdept(ji,jj,jk) - 305) / 460 )  )  &
404                     &              * (-TANH((500 - fsdept(ji,jj,jk)) / 150) + 1) / 2          &
405                     &          + (  35.55 + 1.25 * (5000. - fsdept(ji,jj,jk)) / 5000.         &
406                     &                - 1.62 * TANH( (fsdept(ji,jj,jk) - 60.  ) / 650. )       &
407                     &                + 0.2  * TANH( (fsdept(ji,jj,jk) - 35.  ) / 100. )       &
408                     &                + 0.2  * TANH( (fsdept(ji,jj,jk) - 1000.) / 5000.)    )  &
409                     &              * (-TANH((fsdept(ji,jj,jk) - 500) / 150) + 1) / 2 
410                  tsn(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
411                  tsb(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal)
412               END DO
413            END DO
414         END DO
415         !
416      CASE ( 1 )                  ! T/S data fields read in dta_tem.nc/data_sal.nc files
417         IF(lwp) WRITE(numout,*)
418         IF(lwp) WRITE(numout,*) 'istate_gyre : initial T and S read from dta_tem.nc/data_sal.nc files'
419         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
420         IF(lwp) WRITE(numout,*) '              NetCDF FORMAT'
421
422         ! Read temperature field
423         ! ----------------------
424         CALL iom_open ( 'data_tem', inum )
425         CALL iom_get ( inum, jpdom_data, 'votemper', tsn(:,:,:,jp_tem) ) 
426         CALL iom_close( inum )
427
428         tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:) 
429         tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
430
431         ! Read salinity field
432         ! -------------------
433         CALL iom_open ( 'data_sal', inum )
434         CALL iom_get ( inum, jpdom_data, 'vosaline', tsn(:,:,:,jp_sal) ) 
435         CALL iom_close( inum )
436
437         tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:) 
438         tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
439         !
440      END SELECT
441      !
442      IF(lwp) THEN
443         WRITE(numout,*)
444         WRITE(numout,*) '              Initial temperature and salinity profiles:'
445         WRITE(numout, "(9x,' level   gdept_1d   temperature   salinity   ')" )
446         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 )
447      ENDIF
448      !
449      IF(lwp .AND. lflush) CALL flush(numout)
450      !
451   END SUBROUTINE istate_gyre
452
453
454   SUBROUTINE istate_uvg
455      !!----------------------------------------------------------------------
456      !!                  ***  ROUTINE istate_uvg  ***
457      !!
458      !! ** Purpose :   Compute the geostrophic velocities from (tn,sn) fields
459      !!
460      !! ** Method  :   Using the hydrostatic hypothesis the now hydrostatic
461      !!      pressure is computed by integrating the in-situ density from the
462      !!      surface to the bottom.
463      !!                 p=integral [ rau*g dz ]
464      !!----------------------------------------------------------------------
465      USE dynspg          ! surface pressure gradient             (dyn_spg routine)
466      USE divcur          ! hor. divergence & rel. vorticity      (div_cur routine)
467      USE lbclnk          ! ocean lateral boundary condition (or mpp link)
468      !
469      INTEGER ::   ji, jj, jk        ! dummy loop indices
470      INTEGER ::   indic             ! ???
471      REAL(wp) ::   zmsv, zphv, zmsu, zphu, zalfg     ! temporary scalars
472      REAL(wp), POINTER, DIMENSION(:,:,:) :: zprn
473      !!----------------------------------------------------------------------
474      !
475      CALL wrk_alloc( jpi, jpj, jpk, zprn)
476      !
477      IF(lwp) WRITE(numout,*) 
478      IF(lwp) WRITE(numout,*) 'istate_uvg : Start from Geostrophy'
479      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
480
481      ! Compute the now hydrostatic pressure
482      ! ------------------------------------
483
484      zalfg = 0.5 * grav * rau0
485     
486      zprn(:,:,1) = zalfg * fse3w(:,:,1) * ( 1 + rhd(:,:,1) )       ! Surface value
487
488      DO jk = 2, jpkm1                                              ! Vertical integration from the surface
489         zprn(:,:,jk) = zprn(:,:,jk-1)   &
490            &         + zalfg * fse3w(:,:,jk) * ( 2. + rhd(:,:,jk) + rhd(:,:,jk-1) )
491      END DO 
492
493      ! Compute geostrophic balance
494      ! ---------------------------
495      DO jk = 1, jpkm1
496         DO jj = 2, jpjm1
497            DO ji = fs_2, fs_jpim1   ! vertor opt.
498               zmsv = 1. / MAX(  umask(ji-1,jj+1,jk) + umask(ji  ,jj+1,jk)   &
499                               + umask(ji-1,jj  ,jk) + umask(ji  ,jj  ,jk) , 1.  )
500               zphv = ( zprn(ji  ,jj+1,jk) - zprn(ji-1,jj+1,jk) ) * umask(ji-1,jj+1,jk) / e1u(ji-1,jj+1)   &
501                    + ( zprn(ji+1,jj+1,jk) - zprn(ji  ,jj+1,jk) ) * umask(ji  ,jj+1,jk) / e1u(ji  ,jj+1)   &
502                    + ( zprn(ji  ,jj  ,jk) - zprn(ji-1,jj  ,jk) ) * umask(ji-1,jj  ,jk) / e1u(ji-1,jj  )   &
503                    + ( zprn(ji+1,jj  ,jk) - zprn(ji  ,jj  ,jk) ) * umask(ji  ,jj  ,jk) / e1u(ji  ,jj  )
504               zphv = 1. / rau0 * zphv * zmsv * vmask(ji,jj,jk)
505
506               zmsu = 1. / MAX(  vmask(ji+1,jj  ,jk) + vmask(ji  ,jj  ,jk)   &
507                               + vmask(ji+1,jj-1,jk) + vmask(ji  ,jj-1,jk) , 1.  )
508               zphu = ( zprn(ji+1,jj+1,jk) - zprn(ji+1,jj  ,jk) ) * vmask(ji+1,jj  ,jk) / e2v(ji+1,jj  )   &
509                    + ( zprn(ji  ,jj+1,jk) - zprn(ji  ,jj  ,jk) ) * vmask(ji  ,jj  ,jk) / e2v(ji  ,jj  )   &
510                    + ( zprn(ji+1,jj  ,jk) - zprn(ji+1,jj-1,jk) ) * vmask(ji+1,jj-1,jk) / e2v(ji+1,jj-1)   &
511                    + ( zprn(ji  ,jj  ,jk) - zprn(ji  ,jj-1,jk) ) * vmask(ji  ,jj-1,jk) / e2v(ji  ,jj-1)
512               zphu = 1. / rau0 * zphu * zmsu * umask(ji,jj,jk)
513
514               ! Compute the geostrophic velocities
515               un(ji,jj,jk) = -2. * zphu / ( ff(ji,jj) + ff(ji  ,jj-1) )
516               vn(ji,jj,jk) =  2. * zphv / ( ff(ji,jj) + ff(ji-1,jj  ) )
517            END DO
518         END DO
519      END DO
520
521      IF(lwp) WRITE(numout,*) '         we force to zero bottom velocity'
522
523      ! Susbtract the bottom velocity (level jpk-1 for flat bottom case)
524      ! to have a zero bottom velocity
525
526      DO jk = 1, jpkm1
527         un(:,:,jk) = ( un(:,:,jk) - un(:,:,jpkm1) ) * umask(:,:,jk)
528         vn(:,:,jk) = ( vn(:,:,jk) - vn(:,:,jpkm1) ) * vmask(:,:,jk)
529      END DO
530
531      CALL lbc_lnk( un, 'U', -1. )
532      CALL lbc_lnk( vn, 'V', -1. )
533     
534      ub(:,:,:) = un(:,:,:)
535      vb(:,:,:) = vn(:,:,:)
536     
537      ! WARNING !!!!!
538      ! after initializing u and v, we need to calculate the initial streamfunction bsf.
539      ! Otherwise, only the trend will be computed and the model will blow up (inconsistency).
540      ! to do that, we call dyn_spg with a special trick:
541      ! we fill ua and va with the velocities divided by dt, and the streamfunction will be brought to the
542      ! right value assuming the velocities have been set up in one time step.
543      ! we then set bsfd to zero (first guess for next step is d(psi)/dt = 0.)
544      !  sets up s false trend to calculate the barotropic streamfunction.
545
546      ua(:,:,:) = ub(:,:,:) / rdt
547      va(:,:,:) = vb(:,:,:) / rdt
548
549      ! calls dyn_spg. we assume euler time step, starting from rest.
550      indic = 0
551      CALL dyn_spg( nit000, indic )       ! surface pressure gradient
552
553      ! the new velocity is ua*rdt
554
555      CALL lbc_lnk( ua, 'U', -1. )
556      CALL lbc_lnk( va, 'V', -1. )
557
558      ub(:,:,:) = ua(:,:,:) * rdt
559      vb(:,:,:) = va(:,:,:) * rdt
560      ua(:,:,:) = 0.e0
561      va(:,:,:) = 0.e0
562      un(:,:,:) = ub(:,:,:)
563      vn(:,:,:) = vb(:,:,:)
564       
565      ! Compute the divergence and curl
566
567      CALL div_cur( nit000 )            ! now horizontal divergence and curl
568
569      hdivb(:,:,:) = hdivn(:,:,:)       ! set the before to the now value
570      rotb (:,:,:) = rotn (:,:,:)       ! set the before to the now value
571      !
572      CALL wrk_dealloc( jpi, jpj, jpk, zprn)
573      !
574      IF(lwp .AND. lflush) CALL flush(numout)
575      !
576   END SUBROUTINE istate_uvg
577
578   !!=====================================================================
579END MODULE istate
Note: See TracBrowser for help on using the repository browser.