Changeset 6041 for branches/2015/dev_r5776_UKMO2_OBS_efficiency_improvs/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90

Timestamp:

2015-12-14T10:06:06+01:00 (8 years ago)

Author:

timgraham

Message:

Merged head of trunk into branch

File:

• branches/2015/dev_r5776_UKMO2_OBS_efficiency_improvs/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90 (modified) (13 diffs)

branches/2015/dev_r5776_UKMO2_OBS_efficiency_improvs/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90

@@ -29 +29 @@
    USE daymod          ! calendar
    USE eosbn2          ! eq. of state, Brunt Vaisala frequency (eos     routine)
-   USE ldftra_oce      ! ocean active tracers: lateral physics
+   USE ldftra          ! lateral physics: ocean active tracers
    USE zdf_oce         ! ocean vertical physics
    USE phycst          ! physical constants
    USE dtatsd          ! data temperature and salinity   (dta_tsd routine)
    USE dtauvd          ! data: U & V current             (dta_uvd routine)
-   USE zpshde          ! partial step: hor. derivative (zps_hde routine)
-   USE eosbn2          ! equation of state            (eos bn2 routine)
    USE domvvl          ! varying vertical mesh
-   USE dynspg_oce      ! pressure gradient schemes
-   USE dynspg_flt      ! filtered free surface
-   USE sol_oce         ! ocean solver variables
    !
    USE in_out_manager  ! I/O manager
@@ -76 +71 @@
       !
 
-      IF(lwp) WRITE(numout,*) ' '
+      IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) 'istate_ini : Initialization of the dynamics and tracers'
       IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
 
-      CALL dta_tsd_init                       ! Initialisation of T & S input data
-      IF( lk_c1d ) CALL dta_uvd_init          ! Initialization of U & V input data
+                     CALL dta_tsd_init        ! Initialisation of T & S input data
+      IF( lk_c1d )   CALL dta_uvd_init        ! Initialization of U & V input data
 
       rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk
@@ -103 +98 @@
          ub   (:,:,:) = 0._wp   ;   un   (:,:,:) = 0._wp
          vb   (:,:,:) = 0._wp   ;   vn   (:,:,:) = 0._wp  
-         rotb (:,:,:) = 0._wp   ;   rotn (:,:,:) = 0._wp
-         hdivb(:,:,:) = 0._wp   ;   hdivn(:,:,:) = 0._wp
+                                    hdivn(:,:,:) = 0._wp
          !
          IF( cp_cfg == 'eel' ) THEN
@@ -119 +113 @@
             ENDIF
             IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000
-               CALL wrk_alloc( jpi, jpj, jpk, 2, zuvd )
+               CALL wrk_alloc( jpi,jpj,jpk,2,   zuvd )
                CALL dta_uvd( nit000, zuvd )
                ub(:,:,:) = zuvd(:,:,:,1) ;  un(:,:,:) = ub(:,:,:)
                vb(:,:,:) = zuvd(:,:,:,2) ;  vn(:,:,:) = vb(:,:,:)
-               CALL wrk_dealloc( jpi, jpj, jpk, 2, zuvd )
+               CALL wrk_dealloc( jpi,jpj,jpk,2,   zuvd )
             ENDIF
          ENDIF
-         !
-         CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) )        ! before potential and in situ densities
-#if ! defined key_c1d
-         IF( ln_zps .AND. .NOT. ln_isfcav)                                 &
-            &            CALL zps_hde    ( nit000, jpts, tsb, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
-            &                                            rhd, gru , grv    )  ! of t, s, rd at the last ocean level
-         IF( ln_zps .AND.       ln_isfcav)                                 &
-            &            CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
-            &                                            rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
-            &                                     gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
-#endif
          !   
          ! - ML - sshn could be modified by istate_eel, so that initialization of fse3t_b is done here
@@ -142 +125 @@
             DO jk = 1, jpk
                fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
-            ENDDO
+            END DO
          ENDIF
          ! 
-      ENDIF
-      !
-      IF( lk_agrif ) THEN                  ! read free surface arrays in restart file
-         IF( ln_rstart ) THEN
-            IF( lk_dynspg_flt )  THEN      ! read or initialize the following fields
-               !                           ! gcx, gcxb for agrif_opa_init
-               IF( sol_oce_alloc()  > 0 )   CALL ctl_stop('agrif sol_oce_alloc: allocation of arrays failed')
-               CALL flt_rst( nit000, 'READ' )
-            ENDIF
-         ENDIF                             ! explicit case not coded yet with AGRIF
       ENDIF
       !
@@ -163 +136 @@
       !
       !
-      un_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp
-      ub_b(:,:) = 0._wp ; vb_b(:,:) = 0._wp
+      un_b(:,:) = 0._wp   ;   vn_b(:,:) = 0._wp
+      ub_b(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp
       !
       DO jk = 1, jpkm1
@@ -202 +175 @@
       !! References :  Philander ???
       !!----------------------------------------------------------------------
-      INTEGER  :: ji, jj, jk
-      REAL(wp) ::   zsal = 35.50
+      INTEGER  ::   ji, jj, jk
+      REAL(wp) ::   zsal = 35.50_wp
       !!----------------------------------------------------------------------
       !
@@ -233 +206 @@
       !!                and relative vorticity fields
       !!----------------------------------------------------------------------
-      USE divcur     ! hor. divergence & rel. vorticity      (div_cur routine)
+      USE divhor     ! hor. divergence      (div_hor routine)
       USE iom
       !
@@ -282 +255 @@
             tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
             !
-            ! set the dynamics: U,V, hdiv, rot (and ssh if necessary)
+            ! set the dynamics: U,V, hdiv (and ssh if necessary)
             ! ----------------
             ! Start EEL5 configuration with barotropic geostrophic velocities 
@@ -318 +291 @@
             ENDIF
             !
-            CALL div_cur( nit000 )                  ! horizontal divergence and relative vorticity (curl)
+!!gm  Check  here call to div_hor should not be necessary
+!!gm         div_hor call runoffs  not sure they are defined at that level
+            CALL div_hor( nit000 )                  ! horizontal divergence and relative vorticity (curl)
             ! N.B. the vertical velocity will be computed from the horizontal divergence field
             ! in istate by a call to wzv routine
@@ -371 +346 @@
       !!
       !! ** Method  : - set temprature field
-      !!              - set salinity field
+      !!              - set salinity   field
       !!----------------------------------------------------------------------
       INTEGER :: ji, jj, jk  ! dummy loop indices
@@ -457 +432 @@
       !!                 p=integral [ rau*g dz ]
       !!----------------------------------------------------------------------
-      USE dynspg          ! surface pressure gradient             (dyn_spg routine)
-      USE divcur          ! hor. divergence & rel. vorticity      (div_cur routine)
+      USE divhor          ! hor. divergence                       (div_hor routine)
       USE lbclnk          ! ocean lateral boundary condition (or mpp link)
       !
       INTEGER ::   ji, jj, jk        ! dummy loop indices
-      INTEGER ::   indic             ! ???
       REAL(wp) ::   zmsv, zphv, zmsu, zphu, zalfg     ! temporary scalars
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zprn
       !!----------------------------------------------------------------------
       !
-      CALL wrk_alloc( jpi, jpj, jpk, zprn)
+      CALL wrk_alloc( jpi,jpj,jpk,   zprn)
       !
       IF(lwp) WRITE(numout,*) 
@@ -529 +502 @@
       vb(:,:,:) = vn(:,:,:)
       
-      ! WARNING !!!!!
-      ! after initializing u and v, we need to calculate the initial streamfunction bsf.
-      ! Otherwise, only the trend will be computed and the model will blow up (inconsistency).
-      ! to do that, we call dyn_spg with a special trick:
-      ! we fill ua and va with the velocities divided by dt, and the streamfunction will be brought to the
-      ! right value assuming the velocities have been set up in one time step.
-      ! we then set bsfd to zero (first guess for next step is d(psi)/dt = 0.)
-      !  sets up s false trend to calculate the barotropic streamfunction.
-
-      ua(:,:,:) = ub(:,:,:) / rdt
-      va(:,:,:) = vb(:,:,:) / rdt
-
-      ! calls dyn_spg. we assume euler time step, starting from rest.
-      indic = 0
-      CALL dyn_spg( nit000, indic )       ! surface pressure gradient
-
-      ! the new velocity is ua*rdt
-
-      CALL lbc_lnk( ua, 'U', -1. )
-      CALL lbc_lnk( va, 'V', -1. )
-
-      ub(:,:,:) = ua(:,:,:) * rdt
-      vb(:,:,:) = va(:,:,:) * rdt
-      ua(:,:,:) = 0.e0
-      va(:,:,:) = 0.e0
-      un(:,:,:) = ub(:,:,:)
-      vn(:,:,:) = vb(:,:,:)
-       
-      ! Compute the divergence and curl
-
-      CALL div_cur( nit000 )            ! now horizontal divergence and curl
-
-      hdivb(:,:,:) = hdivn(:,:,:)       ! set the before to the now value
-      rotb (:,:,:) = rotn (:,:,:)       ! set the before to the now value
-      !
-      CALL wrk_dealloc( jpi, jpj, jpk, zprn)
+      !
+!!gm  Check  here call to div_hor should not be necessary
+!!gm         div_hor call runoffs  not sure they are defined at that level
+      CALL div_hor( nit000 )            ! now horizontal divergence
+      !
+      CALL wrk_dealloc( jpi,jpj,jpk,   zprn)
       !
    END SUBROUTINE istate_uvg