Changeset 5111 for branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

Timestamp:

2015-03-02T16:30:57+01:00 (9 years ago)

Author:

mathiot

Message:

add some missing if ln_isfcav, test of option compatibility with ln_isfcav + small documentation

File:

• branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 (modified) (9 diffs)

branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

@@ -298 +298 @@
       ENDIF
 
+!      IF ( ln_isfcav ) THEN
 ! need to be like this to compute the pressure gradient with ISF. If not, level beneath the ISF are not aligned (sum(e3t) /= depth)
 ! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively
-      DO jk = 1, jpkm1
-         e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk) 
-      END DO
-      e3t_1d(jpk) = e3t_1d(jpk-1)   ! we don't care because this level is masked in NEMO
-
-      DO jk = 2, jpk
-         e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1) 
-      END DO
-      e3w_1d(1  ) = 2._wp * (gdept_1d(1) - gdepw_1d(1)) 
+         DO jk = 1, jpkm1
+            e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk) 
+         END DO
+         e3t_1d(jpk) = e3t_1d(jpk-1)   ! we don't care because this level is masked in NEMO
+
+         DO jk = 2, jpk
+            e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1) 
+         END DO
+         e3w_1d(1  ) = 2._wp * (gdept_1d(1) - gdepw_1d(1)) 
+!      END IF
 
 !!gm BUG in s-coordinate this does not work!
@@ -473 +475 @@
          ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code
          IF( cp_cfg == "isomip" .AND. ln_isfcav ) THEN 
-           ! 
-           risfdep(:,:)=200.e0 
-           misfdep(:,:)=1 
-           ij0 = 1 ; ij1 = 40 
-           DO jj = mj0(ij0), mj1(ij1) 
-              risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 
-                END DO 
+            risfdep(:,:)=200.e0 
+            misfdep(:,:)=1 
+            ij0 = 1 ; ij1 = 40 
+            DO jj = mj0(ij0), mj1(ij1) 
+               risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 
+            END DO 
             WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp 
-           ! 
+         ! 
          ELSEIF ( cp_cfg == "isomip2" .AND. ln_isfcav ) THEN
          ! 
@@ -536 +537 @@
             CALL iom_get  ( inum, jpdom_data, 'Bathymetry', bathy )
             CALL iom_close( inum )
-            !  
+            !                                                
             risfdep(:,:)=0._wp         
             misfdep(:,:)=1             
@@ -963 +964 @@
       !!----------------------------------------------------------------------
       !!
-      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
+      INTEGER  ::   ji, jj, jk       ! dummy loop indices
       INTEGER  ::   ik, it           ! temporary integers
-      INTEGER  ::   id, jd, nprocd
       LOGICAL  ::   ll_print         ! Allow  control print for debugging
       REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
       REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t
-      REAL(wp) ::   zmax, zmin       ! Maximum and minimum depth
+      REAL(wp) ::   zmax             ! Maximum depth
       REAL(wp) ::   zdiff            ! temporary scalar
       REAL(wp) ::   zrefdep          ! temporary scalar
@@ -1046 +1046 @@
                      &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   &
                      &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) ))
-                  e3t_0(ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   & 
-                     &                          / ( gdepw_1d(      ik+1) - gdepw_1d(ik) ) 
+                  e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   & 
+                     &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) ) 
                   e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   &
                      &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
@@ -1202 +1202 @@
      
       ! Compute gdep3w_0 (vertical sum of e3w)
-      IF ( ln_isfcav ) THEN
+      IF ( ln_isfcav ) THEN ! if cavity
          WHERE (misfdep == 0) misfdep = 1
          DO jj = 1,jpj
@@ -1335 +1335 @@
             misfdep(jpi,:) = misfdep(  2  ,:) 
          ENDIF
- 
+
          IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN
             mbathy( 1 ,:) = mbathy(jpim1,:)             ! local domain is cyclic east-west
             mbathy(jpi,:) = mbathy(  2  ,:)
          ENDIF
- 
+
          ! split last cell if possible (only where water column is 2 cell or less)
          DO jk = jpkm1, 1, -1
@@ -1358 +1358 @@
             END WHERE
          END DO
- 
+
  
  ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition
@@ -1639 +1639 @@
                IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj+1) ) ibtestjp1 = 0
                ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
-               IF( ibtest == 0 ) THEN
+               IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN
                   mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ;
                END IF