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.
Changeset 4666 for branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Ignore:
Timestamp:
2014-06-11T14:52:23+02:00 (10 years ago)
Author:
mathiot
Message:

#1331 : add ISOMIP config files + ice shelf code

Location:
branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r4488 r4666  
    175175   LOGICAL, PUBLIC ::   ln_zps        !: z-coordinate - partial step 
    176176   LOGICAL, PUBLIC ::   ln_sco        !: s-coordinate or hybrid z-s coordinate 
     177   LOGICAL, PUBLIC ::   ln_isf        !: presence of ISF  
    177178 
    178179   !! All coordinates 
    … …  
    250251   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbkt               !: vertical index of the bottom last T- ocean level 
    251252   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbku, mbkv         !: vertical index of the bottom last U- and W- ocean level 
    252    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bathy              !: ocean depth (meters) 
    253    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_i            !: interior domain T-point mask 
    254    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bmask              !: land/ocean mask of barotropic stream function 
     253   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bathy                              !: ocean depth (meters) 
     254   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_i, umask_i, vmask_i, fmask_i !: interior domain T-point mask 
     255   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bmask                              !: land/ocean mask of barotropic stream function 
     256 
     257   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   micedep                 !: top first ocean level                (ISF) 
     258   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mikt, miku, mikv, mikf  !: first wet T-, U-, V-, F- ocean level (ISF) 
     259   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   icedep, iceh            !: Iceshelf draft, iceshef freeboard    (ISF) 
     260   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   lmask                   !: surface domain T-point mask  
    255261 
    256262   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask   !: land/ocean mask at T-, U-, V- and F-pts 
    … …  
    379385         &      hift  (jpi,jpj) , hifu  (jpi,jpj) , rx1 (jpi,jpj) , STAT=ierr(8) ) 
    380386 
    381       ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) ,                     & 
    382          &     tmask_i(jpi,jpj) , bmask(jpi,jpj) ,                     & 
     387      ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) ,                                      & 
     388         &     tmask_i(jpi,jpj) , umask_i(jpi,jpj), vmask_i(jpi,jpj), fmask_i(jpi,jpj), & 
     389         &     bmask(jpi,jpj)   ,                                                       & 
    383390         &     mbkt   (jpi,jpj) , mbku (jpi,jpj) , mbkv(jpi,jpj) , STAT=ierr(9) ) 
    384391 
     392! (ISF) Allocation of basic array    
     393      ALLOCATE( micedep(jpi,jpj) , icedep(jpi,jpj), iceh(jpi,jpj),      & 
     394         &     mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj) ,           & 
     395         &     mikf(jpi,jpj), lmask(jpi,jpj), STAT=ierr(10) ) 
     396 
    385397      ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk),     &  
    386          &      vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk), STAT=ierr(10) ) 
     398         &      vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk), STAT=ierr(11) ) 
    387399 
    388400#if defined key_noslip_accurate 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r4624 r4666  
    111111      END DO 
    112112      !                                        ! Inverse of the local depth 
    113       hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask(:,:,1) ) * umask(:,:,1) 
    114       hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask(:,:,1) ) * vmask(:,:,1) 
     113      hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:) 
     114      hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:) 
    115115 
    116116                             CALL dom_stp      ! time step 
    … …  
    159159      READ  ( numnam_cfg, namrun, IOSTAT = ios, ERR = 902 ) 
    160160902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp ) 
    161       IF(lwm) WRITE ( numond, namrun ) 
     161      WRITE ( numond, namrun ) 
    162162      ! 
    163163      IF(lwp) THEN                  ! control print 
    … …  
    241241      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 ) 
    242242904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp ) 
    243       IF(lwm) WRITE ( numond, namdom ) 
     243      WRITE ( numond, namdom ) 
    244244 
    245245      IF(lwp) THEN 
    … …  
    303303      READ  ( numnam_cfg, namcla, IOSTAT = ios, ERR = 906 ) 
    304304906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in configuration namelist', lwp ) 
    305       IF(lwm) WRITE( numond, namcla ) 
     305      WRITE( numond, namcla ) 
    306306 
    307307      IF(lwp) THEN 
    … …  
    327327      READ  ( numnam_cfg, namnc4, IOSTAT = ios, ERR = 908 ) 
    328328908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in configuration namelist', lwp ) 
    329       IF(lwm) WRITE( numond, namnc4 ) 
     329      WRITE( numond, namnc4 ) 
    330330 
    331331      IF(lwp) THEN                        ! control print 
    … …  
    365365      ! 
    366366      IF(lk_mpp) THEN 
    367          CALL mpp_minloc( e1t(:,:), tmask(:,:,1), ze1min, iimi1,ijmi1 ) 
    368          CALL mpp_minloc( e2t(:,:), tmask(:,:,1), ze2min, iimi2,ijmi2 ) 
    369          CALL mpp_maxloc( e1t(:,:), tmask(:,:,1), ze1max, iima1,ijma1 ) 
    370          CALL mpp_maxloc( e2t(:,:), tmask(:,:,1), ze2max, iima2,ijma2 ) 
     367         CALL mpp_minloc( e1t(:,:), tmask_i(:,:), ze1min, iimi1,ijmi1 ) 
     368         CALL mpp_minloc( e2t(:,:), tmask_i(:,:), ze2min, iimi2,ijmi2 ) 
     369         CALL mpp_maxloc( e1t(:,:), tmask_i(:,:), ze1max, iima1,ijma1 ) 
     370         CALL mpp_maxloc( e2t(:,:), tmask_i(:,:), ze2max, iima2,ijma2 ) 
    371371      ELSE 
    372          ze1min = MINVAL( e1t(:,:), mask = tmask(:,:,1) == 1._wp )     
    373          ze2min = MINVAL( e2t(:,:), mask = tmask(:,:,1) == 1._wp )     
    374          ze1max = MAXVAL( e1t(:,:), mask = tmask(:,:,1) == 1._wp )     
    375          ze2max = MAXVAL( e2t(:,:), mask = tmask(:,:,1) == 1._wp )     
    376  
    377          iloc  = MINLOC( e1t(:,:), mask = tmask(:,:,1) == 1._wp ) 
     372         ze1min = MINVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )     
     373         ze2min = MINVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )     
     374         ze1max = MAXVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )     
     375         ze2max = MAXVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )     
     376 
     377         iloc  = MINLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp ) 
    378378         iimi1 = iloc(1) + nimpp - 1 
    379379         ijmi1 = iloc(2) + njmpp - 1 
    380          iloc  = MINLOC( e2t(:,:), mask = tmask(:,:,1) == 1._wp ) 
     380         iloc  = MINLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp ) 
    381381         iimi2 = iloc(1) + nimpp - 1 
    382382         ijmi2 = iloc(2) + njmpp - 1 
    383          iloc  = MAXLOC( e1t(:,:), mask = tmask(:,:,1) == 1._wp ) 
     383         iloc  = MAXLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp ) 
    384384         iima1 = iloc(1) + nimpp - 1 
    385385         ijma1 = iloc(2) + njmpp - 1 
    386          iloc  = MAXLOC( e2t(:,:), mask = tmask(:,:,1) == 1._wp ) 
     386         iloc  = MAXLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp ) 
    387387         iima2 = iloc(1) + nimpp - 1 
    388388         ijma2 = iloc(2) + njmpp - 1 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90

    r4366 r4666  
    638638      CALL iom_close( inum ) 
    639639       
     640! need to be define for the extended grid south of -80S 
     641! some point are undefined but you need to have e1 and e2 .NE. 0 
     642      WHERE (e1t==0.0_wp) 
     643         e1t=1.0e2 
     644      END WHERE 
     645      WHERE (e1v==0.0_wp) 
     646         e1v=1.0e2 
     647      END WHERE 
     648      WHERE (e1u==0.0_wp) 
     649         e1u=1.0e2 
     650      END WHERE 
     651      WHERE (e1f==0.0_wp) 
     652         e1f=1.0e2 
     653      END WHERE 
     654      WHERE (e2t==0.0_wp) 
     655         e2t=1.0e2 
     656      END WHERE 
     657      WHERE (e2v==0.0_wp) 
     658         e2v=1.0e2 
     659      END WHERE 
     660      WHERE (e2u==0.0_wp) 
     661         e2u=1.0e2 
     662      END WHERE 
     663      WHERE (e2f==0.0_wp) 
     664         e2f=1.0e2 
     665      END WHERE 
     666        
    640667    END SUBROUTINE hgr_read 
    641668     

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90

    r4624 r4666  
    184184         END DO   
    185185      END DO   
     186       
     187      ! (ISF) define barotropic mask and mask the ice shelf point 
     188      lmask(:,:)=tmask(:,:,1) ! at this stage ice shelf is not masked 
     189       
     190      DO jk = 1, jpk 
     191         DO jj = 1, jpj 
     192            DO ji = 1, jpi 
     193               IF( REAL( micedep(ji,jj) - jk, wp ) - 0.1_wp >= 0._wp )   THEN 
     194                  tmask(ji,jj,jk) = 0._wp 
     195               END IF 
     196            END DO   
     197         END DO   
     198      END DO   
    186199 
    187200!!gm  ???? 
    … …  
    207220      ! Interior domain mask (used for global sum) 
    208221      ! -------------------- 
    209       tmask_i(:,:) = tmask(:,:,1) 
     222      tmask_i(:,:) = lmask(:,:)            ! (ISH) tmask_i = 1 even on the ice shelf 
    210223      iif = jpreci                         ! ??? 
    211224      iil = nlci - jpreci + 1 
    … …  
    250263         END DO 
    251264      END DO 
     265      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet u point 
     266      DO jj = 1, jpjm1 
     267         DO ji = 1, fs_jpim1   ! vector loop 
     268            umask_i(ji,jj)  = lmask(ji,jj) * lmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:))) 
     269            vmask_i(ji,jj)  = lmask(ji,jj) * lmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:))) 
     270         END DO 
     271         DO ji = 1, jpim1      ! NO vector opt. 
     272            fmask_i(ji,jj) =  lmask(ji,jj  ) * lmask(ji+1,jj  )   & 
     273               &            * lmask(ji,jj+1) * lmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:))) 
     274         END DO 
     275      END DO 
    252276      CALL lbc_lnk( umask, 'U', 1._wp )      ! Lateral boundary conditions 
    253277      CALL lbc_lnk( vmask, 'V', 1._wp ) 
    254278      CALL lbc_lnk( fmask, 'F', 1._wp ) 
     279      CALL lbc_lnk( umask_i, 'U', 1._wp )      ! Lateral boundary conditions 
     280      CALL lbc_lnk( vmask_i, 'V', 1._wp ) 
     281      CALL lbc_lnk( fmask_i, 'F', 1._wp ) 
    255282 
    256283 
    257284      ! 4. ocean/land mask for the elliptic equation 
    258285      ! -------------------------------------------- 
    259       bmask(:,:) = tmask(:,:,1)       ! elliptic equation is written at t-point 
     286      bmask(:,:) = lmask(:,:)       ! elliptic equation is written at t-point 
    260287      ! 
    261288      !                               ! Boundary conditions 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r4624 r4666  
    169169      fsdept_b(:,:,1) = 0.5_wp * fse3w_b(:,:,1) 
    170170      fsdepw_b(:,:,1) = 0.0_wp 
    171       DO jk = 2, jpk 
    172          fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk) 
    173          fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1) 
    174          fsde3w_n(:,:,jk) = fsdept_n(:,:,jk  ) - sshn   (:,:) 
    175          fsdept_b(:,:,jk) = fsdept_b(:,:,jk-1) + fse3w_b(:,:,jk) 
    176          fsdepw_b(:,:,jk) = fsdepw_b(:,:,jk-1) + fse3t_b(:,:,jk-1) 
     171      DO jj = 1,jpj 
     172         DO ji = 1,jpi 
     173            DO jk = 1,mikt(ji,jj)-1 
     174               fsdept_n(ji,jj,jk) = gdept_0(ji,jj,jk) 
     175               fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 
     176               fsde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj) 
     177               fsdept_b(ji,jj,jk) = gdept_0(ji,jj,jk) 
     178               fsdepw_b(ji,jj,jk) = gdepw_0(ji,jj,jk) 
     179            END DO 
     180            DO jk = mikt(ji,jj), jpk 
     181               fsdept_n(ji,jj,jk) = fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk) 
     182               fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 
     183               fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj) 
     184               fsdept_b(ji,jj,jk) = fsdept_b(ji,jj,jk-1) + fse3w_b(ji,jj,jk) 
     185               fsdepw_b(ji,jj,jk) = fsdepw_b(ji,jj,jk-1) + fse3t_b(ji,jj,jk-1) 
     186            END DO 
     187         END DO 
    177188      END DO 
    178189 
    … …  
    185196         hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 
    186197      END DO 
    187       hur_b(:,:) = umask(:,:,1) / ( hu_b(:,:) + 1. - umask(:,:,1) ) 
    188       hvr_b(:,:) = vmask(:,:,1) / ( hv_b(:,:) + 1. - vmask(:,:,1) ) 
     198      hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1. - umask_i(:,:) ) 
     199      hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1. - vmask_i(:,:) ) 
    189200 
    190201      ! Restoring frequencies for z_tilde coordinate 
    … …  
    293304      !                                                ! --------------------------------------------- ! 
    294305 
    295       z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask(:,:,1) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) ) 
     306      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask_i(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask_i(:,:) ) 
    296307      DO jk = 1, jpkm1 
    297308         ! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0) 
    … …  
    313324            zht  (:,:) = zht  (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 
    314325         END DO 
    315          zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask(:,:,1) ) 
     326         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) 
    316327 
    317328         ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only) 
    … …  
    338349         ELSE                         ! layer case 
    339350            DO jk = 1, jpkm1 
    340                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) 
     351               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 
    341352            END DO 
    342353         END IF 
    … …  
    463474      !                                           ! ---baroclinic part--------- ! 
    464475         DO jk = 1, jpkm1 
    465             fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) 
     476            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 
    466477         END DO 
    467478      ENDIF 
    … …  
    531542      END DO 
    532543      !                                        ! Inverse of the local depth 
    533       hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask(:,:,1) ) * umask(:,:,1) 
    534       hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask(:,:,1) ) * vmask(:,:,1) 
     544      hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:) 
     545      hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:) 
    535546 
    536547      CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) 
    … …  
    830841               DO jk=1,jpk 
    831842                  fse3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & 
    832                       &                            / ( ht_0(:,:) + 1._wp - tmask(:,:,1) ) * tmask(:,:,jk) & 
     843                      &                            / ( ht_0(:,:) + 1._wp - tmask_i(:,:) ) * tmask(:,:,jk) & 
    833844                      &            + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk)) 
    834845               END DO 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90

    r4292 r4666  
    132132       
    133133      CALL dom_uniq( zprw, 'T' ) 
    134       zprt = tmask(:,:,1) * zprw                               !    ! unique point mask 
     134      DO jj = 1, jpj 
     135         DO ji = 1, jpi 
     136            jk=mikt(ji,jj)  
     137            zprt(ji,jj) = tmask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask 
     138         END DO 
     139      END DO                             !    ! unique point mask 
    135140      CALL iom_rstput( 0, 0, inum2, 'tmaskutil', zprt, ktype = jp_i1 )   
    136141      CALL dom_uniq( zprw, 'U' ) 
    137       zprt = umask(:,:,1) * zprw 
     142      DO jj = 1, jpj 
     143         DO ji = 1, jpi 
     144            jk=miku(ji,jj)  
     145            zprt(ji,jj) = umask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask 
     146         END DO 
     147      END DO 
    138148      CALL iom_rstput( 0, 0, inum2, 'umaskutil', zprt, ktype = jp_i1 )   
    139149      CALL dom_uniq( zprw, 'V' ) 
    140       zprt = vmask(:,:,1) * zprw 
     150      DO jj = 1, jpj 
     151         DO ji = 1, jpi 
     152            jk=mikv(ji,jj)  
     153            zprt(ji,jj) = vmask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask 
     154         END DO 
     155      END DO 
    141156      CALL iom_rstput( 0, 0, inum2, 'vmaskutil', zprt, ktype = jp_i1 )   
    142157      CALL dom_uniq( zprw, 'F' ) 
    143       zprt = fmask(:,:,1) * zprw 
     158      DO jj = 1, jpj 
     159         DO ji = 1, jpi 
     160            jk=mikf(ji,jj)  
     161            zprt(ji,jj) = fmask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask 
     162         END DO 
     163      END DO 
    144164      CALL iom_rstput( 0, 0, inum2, 'fmaskutil', zprt, ktype = jp_i1 )   
    145165 
    … …  
    168188       
    169189      ! note that mbkt is set to 1 over land ==> use surface tmask 
    170       zprt(:,:) = tmask(:,:,1) * REAL( mbkt(:,:) , wp ) 
     190      zprt(:,:) = lmask(:,:) * REAL( mbkt(:,:) , wp ) 
    171191      CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i2 )     !    ! nb of ocean T-points 
     192      zprt(:,:) = lmask(:,:) * REAL( mikt(:,:) , wp ) 
     193      CALL iom_rstput( 0, 0, inum4, 'misf', zprt, ktype = jp_i2 )       !    ! nb of ocean T-points 
     194      zprt(:,:) = lmask(:,:) * REAL( icedep(:,:) , wp ) 
     195      CALL iom_rstput( 0, 0, inum4, 'isfdraft', zprt, ktype = jp_r4 )       !    ! nb of ocean T-points 
    172196             
    173197      IF( ln_sco ) THEN                                         ! s-coordinate 
    … …  
    203227            DO jj = 1,jpj    
    204228               DO ji = 1,jpi 
    205                   e3tp(ji,jj) = e3t_0(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1) 
    206                   e3wp(ji,jj) = e3w_0(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1) 
     229                  e3tp(ji,jj) = e3t_0(ji,jj,mbkt(ji,jj)) * lmask(ji,jj) 
     230                  e3wp(ji,jj) = e3w_0(ji,jj,mbkt(ji,jj)) * lmask(ji,jj) 
    207231               END DO 
    208232            END DO 
    … …  
    228252            DO jj = 1,jpj    
    229253               DO ji = 1,jpi 
    230                   zprt(ji,jj) = gdept_0(ji,jj,mbkt(ji,jj)  ) * tmask(ji,jj,1) 
    231                   zprw(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) * tmask(ji,jj,1) 
     254                  zprt(ji,jj) = gdept_0(ji,jj,mbkt(ji,jj)  ) * lmask(ji,jj) 
     255                  zprw(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) * lmask(ji,jj) 
    232256               END DO 
    233257            END DO 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r4624 r4666  
    3535   USE oce               ! ocean variables 
    3636   USE dom_oce           ! ocean domain 
     37   USE sbc_oce           ! surface variable (isf) 
    3738   USE closea            ! closed seas 
    3839   USE c1d               ! 1D vertical configuration 
    … …  
    102103      ! 
    103104      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco 
     105      NAMELIST/namsbc/ nn_fsbc   , ln_ana , ln_flx  , ln_blk_clio, ln_blk_core, ln_cpl,   & 
     106         &             ln_blk_mfs, ln_apr_dyn, nn_ice , ln_dm2dc, ln_rnf, ln_ssr, nn_fwb, ln_cdgw, nn_isf 
    104107      !!---------------------------------------------------------------------- 
    105108      ! 
    … …  
    145148      IF( .NOT.lk_c1d )   CALL zgr_bat_ctl      ! check bathymetry (mbathy) and suppress isolated ocean points 
    146149                          CALL zgr_bot_level    ! deepest ocean level for t-, u- and v-points 
     150                          CALL zgr_top_level    ! shallowest ocean level for T-, U-, V- points 
    147151      ! 
    148152      IF( lk_c1d ) THEN                         ! 1D config.: same mbathy value over the 3x3 domain 
    … …  
    294298         gdepw_1d(1) = 0._wp                    ! force first w-level to be exactly at zero 
    295299      ENDIF 
     300 
     301! need to be like this to compute the pressure gradient with ISF 
     302! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively 
     303      DO jk = 1, jpkm1 
     304         e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk)  
     305      END DO 
     306      e3t_1d(jpk) = e3t_1d(jpk-1)   ! we don't care because this level is masked in NEMO 
     307 
     308      DO jk = 2, jpk 
     309         e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1)  
     310      END DO 
     311      e3w_1d(1  ) = 2._wp * (gdept_1d(1) - gdepw_1d(1))  
    296312 
    297313!!gm BUG in s-coordinate this does not work! 
    … …  
    451467         END DO 
    452468         ! 
     469         ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code 
     470         IF( cp_cfg == "isomip" ) THEN  
     471           !  
     472           icedep(:,:)=200.e0  
     473           micedep(:,:)=1  
     474           ij0 = 1 ; ij1 = 40  
     475           DO jj = mj0(ij0), mj1(ij1)  
     476              icedep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp  
     477                END DO  
     478            WHERE( bathy(:,:) <= 0._wp )  icedep(:,:) = 0._wp  
     479           !  
     480         ELSEIF ( cp_cfg == "isomip2" ) THEN 
     481         !  
     482            icedep(:,:)=0.e0 
     483            micedep(:,:)=1 
     484            ij0 = 1 ; ij1 = 40 
     485            DO jj = mj0(ij0), mj1(ij1) 
     486               icedep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 
     487            END DO 
     488            WHERE( bathy(:,:) <= 0._wp )  icedep(:,:) = 0._wp 
     489         END IF 
     490         ! 
    453491         !                                            ! ================ ! 
    454492      ELSEIF( ntopo == 1 ) THEN                       !   read in file   ! (over the local domain) 
    … …  
    492530            CALL iom_get  ( inum, jpdom_data, 'Bathymetry', bathy ) 
    493531            CALL iom_close( inum ) 
    494             !                                                 
     532            !   
     533            icedep(:,:)=0._wp          
     534            micedep(:,:)=1              
     535            IF ( nn_isf == 1 .OR. nn_isf == 4 ) THEN 
     536               CALL iom_open ( 'isf_draft_meter.nc', inum )  
     537               CALL iom_get  ( inum, jpdom_data, 'isf_draft', icedep ) 
     538               CALL iom_close( inum ) 
     539               WHERE( bathy(:,:) <= 0._wp )  icedep(:,:) = 0._wp 
     540            END IF 
     541            !        
    495542            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration 
    496543               ! 
    … …  
    783830   END SUBROUTINE zgr_bot_level 
    784831 
     832      SUBROUTINE zgr_top_level 
     833      !!---------------------------------------------------------------------- 
     834      !!                    ***  ROUTINE zgr_bot_level  *** 
     835      !! 
     836      !! ** Purpose :   defines the vertical index of ocean top (mik. arrays) 
     837      !! 
     838      !! ** Method  :   computes from micedep with a minimum value of 1 
     839      !! 
     840      !! ** Action  :   mikt, miku, mikv :   vertical indices of the shallowest  
     841      !!                                     ocean level at t-, u- & v-points 
     842      !!                                     (min value = 1) 
     843      !!---------------------------------------------------------------------- 
     844      !! 
     845      INTEGER ::   ji, jj   ! dummy loop indices 
     846      REAL(wp), POINTER, DIMENSION(:,:) ::  zmik 
     847      !!---------------------------------------------------------------------- 
     848      ! 
     849      IF( nn_timing == 1 )  CALL timing_start('zgr_top_level') 
     850      ! 
     851      CALL wrk_alloc( jpi, jpj, zmik ) 
     852      ! 
     853      IF(lwp) WRITE(numout,*) 
     854      IF(lwp) WRITE(numout,*) '    zgr_top_level : ocean top k-index of T-, U-, V- and W-levels ' 
     855      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~' 
     856      ! 
     857      mikt(:,:) = MAX( micedep(:,:) , 1 )    ! top k-index of T-level (=1) 
     858      !                                      ! top k-index of W-level (=mikt) 
     859      DO jj = 1, jpjm1                       ! top k-index of U- (U-) level 
     860         DO ji = 1, jpim1 
     861            miku(ji,jj) = MAX(  mikt(ji+1,jj  ) , mikt(ji,jj)  ) 
     862            mikv(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj)  ) 
     863            mikf(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj), mikt(ji+1,jj  ), mikt(ji+1,jj+1)  ) 
     864         END DO 
     865      END DO 
     866 
     867      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk  
     868      zmik(:,:) = REAL( miku(:,:), wp )   ;   CALL lbc_lnk(zmik,'U',1.)   ;   miku  (:,:) = MAX( INT( zmik(:,:) ), 1 ) 
     869      zmik(:,:) = REAL( mikv(:,:), wp )   ;   CALL lbc_lnk(zmik,'V',1.)   ;   mikv  (:,:) = MAX( INT( zmik(:,:) ), 1 ) 
     870      zmik(:,:) = REAL( mikf(:,:), wp )   ;   CALL lbc_lnk(zmik,'F',1.)   ;   mikf  (:,:) = MAX( INT( zmik(:,:) ), 1 ) 
     871      ! 
     872      CALL wrk_dealloc( jpi, jpj, zmik ) 
     873      ! 
     874      IF( nn_timing == 1 )  CALL timing_stop('zgr_top_level') 
     875      ! 
     876   END SUBROUTINE zgr_top_level 
    785877 
    786878   SUBROUTINE zgr_zco 
    … …  
    861953      !!---------------------------------------------------------------------- 
    862954      !! 
    863       INTEGER  ::   ji, jj, jk       ! dummy loop indices 
     955      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    864956      INTEGER  ::   ik, it           ! temporary integers 
     957      INTEGER  ::   id, jd, nprocd 
     958      INTEGER  ::   icompt, ibtest   ! (ISF) 
    865959      LOGICAL  ::   ll_print         ! Allow  control print for debugging 
    866960      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points 
    867961      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t 
    868       REAL(wp) ::   zmax             ! Maximum depth 
     962      REAL(wp) ::   zmax, zmin       ! Maximum and minimum depth 
    869963      REAL(wp) ::   zdiff            ! temporary scalar 
    870964      REAL(wp) ::   zrefdep          ! temporary scalar 
     965      REAL(wp) ::   eps=0.99            ! small offset to avoid large pool in case bathy slightly greater than icedep 
     966      REAL(wp), POINTER, DIMENSION(:,:)   ::   zbathy, zmask   ! 3D workspace (ISH) 
     967      INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmicedep   ! 3D workspace (ISH) 
    871968      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zprt 
    872969      !!--------------------------------------------------------------------- 
    … …  
    875972      ! 
    876973      CALL wrk_alloc( jpi, jpj, jpk, zprt ) 
     974      CALL wrk_alloc( jpi, jpj, zbathy, zmask) 
     975      CALL wrk_alloc( jpi, jpj, zmbathy, zmicedep) 
    877976      ! 
    878977      IF(lwp) WRITE(numout,*) 
    … …  
    9061005         WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1 
    9071006      END DO 
     1007      ! (ISF) compute micedep 
     1008      WHERE( icedep(:,:) == 0._wp ) ;   micedep(:,:) = 1   ! no ice shelf : set micedep to 1   
     1009      ELSEWHERE                     ;   micedep(:,:) = 2   ! iceshelf : initialize micedep to second level  
     1010      END WHERE   
     1011 
     1012      ! Compute micedep for ocean points (i.e. first wet level)  
     1013      ! find the first ocean level such that the first level thickness  
     1014      ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where  
     1015      ! e3t_0 is the reference level thickness  
     1016      DO jk = 2, jpkm1  
     1017         zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )  
     1018         WHERE( 0._wp < icedep(:,:) .AND. icedep(:,:) >= zdepth )   micedep(:,:) = jk+1  
     1019      END DO  
     1020      WHERE (icedep <= e3t_1d(1) .AND. icedep .GT. 0._wp) 
     1021         icedep = 0. 
     1022         micedep= 1 
     1023      END WHERE 
     1024     
     1025 
     1026      icompt = 0  
     1027      DO jl = 1, 10  
     1028        IF( lk_mpp ) THEN 
     1029           zbathy(:,:) = FLOAT( micedep(:,:) ) 
     1030           CALL lbc_lnk( zbathy, 'T', 1. ) 
     1031           micedep(:,:) = INT( zbathy(:,:) ) 
     1032           CALL lbc_lnk( icedep, 'T', 1. ) 
     1033           CALL lbc_lnk( bathy, 'T', 1. ) 
     1034           zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1035           CALL lbc_lnk( zbathy, 'T', 1. ) 
     1036           mbathy(:,:) = INT( zbathy(:,:) ) 
     1037        ENDIF 
     1038         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN  
     1039            micedep( 1 ,:) = micedep(jpim1,:)           ! local domain is cyclic east-west  
     1040            micedep(jpi,:) = micedep(  2  ,:)  
     1041         ENDIF 
     1042 
     1043         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN 
     1044            mbathy( 1 ,:) = mbathy(jpim1,:)             ! local domain is cyclic east-west 
     1045            mbathy(jpi,:) = mbathy(  2  ,:) 
     1046         ENDIF 
     1047 
     1048         WHERE (mbathy == 0)  
     1049            icedep = 0._wp 
     1050            micedep= 0 
     1051            bathy  = 0._wp 
     1052         ENDWHERE 
     1053 
     1054         WHERE (bathy(:,:) < icedep(:,:)+eps) 
     1055            micedep(:,:) = 0  
     1056            icedep(:,:)  = 0._wp 
     1057            mbathy(:,:)  = 0 
     1058            bathy(:,:)   = 0._wp 
     1059         END WHERE 
     1060 
     1061         DO jj = 1, jpj 
     1062            DO ji = 1, jpi 
     1063               IF (bathy(ji,jj) .GT. icedep(ji,jj) .AND. mbathy(ji,jj) .LT. micedep(ji,jj)) THEN 
     1064                  bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 
     1065                  mbathy(ji,jj)= mbathy(ji,jj) + 1 
     1066               END IF 
     1067            END DO 
     1068         END DO 
     1069 
     1070 
     1071         ! At least 2 levels for water thickness at T, U, and V point. 
     1072        zmicedep(:,:)=micedep(:,:) 
     1073        zmbathy (:,:)=mbathy (:,:) 
     1074 
     1075         DO jj = 2, jpjm1 
     1076            DO ji = 2, jpim1 
     1077               IF( zmicedep(ji,jj) == zmbathy(ji,jj) .AND. zmbathy(ji,jj) .GT. 1) THEN 
     1078                  mbathy(ji,jj) = zmbathy(ji,jj) + 1  
     1079                  bathy(ji,jj)=gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj))*e3zps_rat ) 
     1080               ENDIF 
     1081               IF( zmicedep(ji,jj+1) == zmbathy(ji,jj) .AND. zmbathy(ji,jj) .GT. 1) THEN 
     1082                  mbathy(ji,jj) = zmbathy(ji,jj) + 1  
     1083                  bathy(ji,jj)=gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj))*e3zps_rat ) 
     1084               ENDIF 
     1085               IF( zmicedep(ji,jj-1) == zmbathy(ji,jj) .AND. zmbathy(ji,jj) .GT. 1) THEN 
     1086                  mbathy(ji,jj) = zmbathy(ji,jj) + 1  
     1087                  bathy(ji,jj)=gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj))*e3zps_rat ) 
     1088               ENDIF 
     1089               IF( zmicedep(ji+1,jj) == zmbathy(ji,jj) .AND. zmbathy(ji,jj) .GT. 1) THEN 
     1090                  mbathy(ji,jj) = zmbathy(ji,jj) + 1  
     1091                  bathy(ji,jj)=gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj))*e3zps_rat ) 
     1092               ENDIF 
     1093               IF( zmicedep(ji-1,jj) == zmbathy(ji,jj) .AND. zmbathy(ji,jj) .GT. 1) THEN 
     1094                  mbathy(ji,jj) = zmbathy(ji,jj) + 1  
     1095                  bathy(ji,jj)=gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj))*e3zps_rat ) 
     1096               ENDIF 
     1097            ENDDO 
     1098         ENDDO 
     1099   IF( lk_mpp ) THEN  
     1100      zbathy(:,:) = FLOAT( micedep(:,:) )  
     1101      CALL lbc_lnk( zbathy, 'T', 1. )  
     1102      micedep(:,:) = INT( zbathy(:,:) )  
     1103      CALL lbc_lnk( icedep, 'T', 1. )  
     1104           CALL lbc_lnk( bathy, 'T', 1. ) 
     1105           zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1106           CALL lbc_lnk( zbathy, 'T', 1. ) 
     1107           mbathy(:,:) = INT( zbathy(:,:) ) 
     1108   ENDIF  
     1109 
     1110         ! if single ocean point put as land 
     1111         zmask=1 
     1112         WHERE (mbathy .EQ. 0) zmask=0 
     1113         DO jj = 2, jpjm1 
     1114            DO ji = 2, jpim1 
     1115               ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 
     1116               IF (ibtest .LE. 1) THEN 
     1117                  bathy(ji,jj)=0._wp 
     1118                  mbathy(ji,jj)=0 
     1119                  icedep(ji,jj)=0._wp 
     1120                  micedep(ji,jj)=0 
     1121               END IF 
     1122            END DO 
     1123         END DO 
     1124   IF( lk_mpp ) THEN  
     1125      zbathy(:,:) = FLOAT( micedep(:,:) )  
     1126      CALL lbc_lnk( zbathy, 'T', 1. )  
     1127      micedep(:,:) = INT( zbathy(:,:) )  
     1128      CALL lbc_lnk( icedep, 'T', 1. )  
     1129           CALL lbc_lnk( bathy, 'T', 1. ) 
     1130           zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1131           CALL lbc_lnk( zbathy, 'T', 1. ) 
     1132           mbathy(:,:) = INT( zbathy(:,:) ) 
     1133   ENDIF  
     1134 
     1135         ! if single point on isf coast line 
     1136         DO jk = 1, jpk 
     1137         WHERE (micedep==0) micedep=jpk 
     1138         zmask=0 
     1139         WHERE (micedep .LE. jk) zmask=1 
     1140         DO jj = 2, jpjm1 
     1141            DO ji = 2, jpim1 
     1142               IF (micedep(ji,jj) .EQ. jk) THEN 
     1143               ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 
     1144               IF (ibtest .LE. 1) THEN 
     1145                     icedep(ji,jj)=gdepw_1d(jk+1) ; micedep(ji,jj)=jk+1 
     1146                     IF (micedep(ji,jj) .GT. mbathy(ji,jj)) micedep(ji,jj) = jpk 
     1147                !     bathy(ji,jj)=0.  ; mbathy(ji,jj)=0 
     1148                  !END IF 
     1149               END IF 
     1150               END IF 
     1151            END DO 
     1152         END DO 
     1153         WHERE (micedep==jpk)  
     1154             micedep=0 ; icedep=0._wp ; mbathy=0 ; bathy=0._wp 
     1155         END WHERE 
     1156 
     1157   IF( lk_mpp ) THEN  
     1158      zbathy(:,:) = FLOAT( micedep(:,:) )  
     1159      CALL lbc_lnk( zbathy, 'T', 1. )  
     1160      micedep(:,:) = INT( zbathy(:,:) )  
     1161      CALL lbc_lnk( icedep, 'T', 1. )  
     1162           CALL lbc_lnk( bathy, 'T', 1. ) 
     1163           zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1164           CALL lbc_lnk( zbathy, 'T', 1. ) 
     1165           mbathy(:,:) = INT( zbathy(:,:) ) 
     1166   ENDIF  
     1167         END DO 
     1168 
     1169! fill hole in ice shelf 
     1170         WHERE (micedep==0) micedep=jpk 
     1171        zmicedep(:,:)=micedep(:,:) 
     1172        zmbathy (:,:)=mbathy (:,:) 
     1173         DO jj = 2, jpjm1 
     1174            DO ji = 2, jpim1 
     1175                  ibtest=MIN(zmicedep(ji-1,jj), zmicedep(ji+1,jj), zmicedep(ji,jj-1), zmicedep(ji,jj+1)) 
     1176                  IF( ibtest > zmicedep(ji,jj)) THEN 
     1177                     micedep(ji,jj) = ibtest 
     1178                     icedep(ji,jj)  = gdepw_1d(ibtest) 
     1179                  ENDIF 
     1180            ENDDO 
     1181         ENDDO 
     1182         WHERE (micedep==jpk)  
     1183             micedep=0 ; icedep=0. ; mbathy=0 ; bathy=0 
     1184         END WHERE 
     1185         IF( lk_mpp ) THEN  
     1186            zbathy(:,:) = FLOAT( micedep(:,:) )  
     1187            CALL lbc_lnk( zbathy, 'T', 1. )  
     1188            micedep(:,:) = INT( zbathy(:,:) )  
     1189            CALL lbc_lnk( icedep, 'T', 1. )  
     1190            CALL lbc_lnk( bathy, 'T', 1. ) 
     1191            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1192            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1193            mbathy(:,:) = INT( zbathy(:,:) ) 
     1194        ENDIF  
     1195 
     1196! fill hole in bathymetry 
     1197        zmicedep(:,:)=micedep(:,:) 
     1198        zmbathy (:,:)=mbathy (:,:) 
     1199         DO jj = 2, jpjm1 
     1200            DO ji = 2, jpim1 
     1201               ibtest = MAX(  zmbathy(ji-1,jj), zmbathy(ji+1,jj),   & 
     1202                  &           zmbathy(ji,jj-1), zmbathy(ji,jj+1)  ) 
     1203               IF( ibtest < zmbathy(ji,jj) ) THEN 
     1204                  mbathy(ji,jj) = ibtest 
     1205                  bathy(ji,jj)  = gdepw_1d(ibtest+1)  
     1206               ENDIF 
     1207            END DO 
     1208         END DO 
     1209         IF( lk_mpp ) THEN  
     1210            zbathy(:,:) = FLOAT( micedep(:,:) )  
     1211            CALL lbc_lnk( zbathy, 'T', 1. )  
     1212            micedep(:,:) = INT( zbathy(:,:) )  
     1213            CALL lbc_lnk( icedep, 'T', 1. )  
     1214            CALL lbc_lnk( bathy, 'T', 1. ) 
     1215            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1216            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1217            mbathy(:,:) = INT( zbathy(:,:) ) 
     1218        ENDIF  
     1219        ! remove pool of water stuck between ice shelf and bathymetry 
     1220        DO jk = 1, jpk 
     1221        WHERE (micedep==0) micedep=jpk 
     1222        zmicedep(:,:)=micedep(:,:) 
     1223        zmbathy (:,:)=mbathy (:,:) 
     1224        DO jj = 2, jpjm1 
     1225            DO ji = 2, jpim1 
     1226               IF( jk .GE. zmicedep(ji,jj) .AND. jk .LE. zmbathy(ji,jj) ) THEN 
     1227               IF( (jk > zmbathy(ji,jj+1) .OR. jk < zmicedep(ji,jj+1)) .AND. & 
     1228                 & (jk > zmbathy(ji,jj-1) .OR. jk < zmicedep(ji,jj-1)) .AND. & 
     1229                 & (jk > zmbathy(ji+1,jj) .OR. jk < zmicedep(ji+1,jj)) .AND. & 
     1230                 & (jk > zmbathy(ji-1,jj) .OR. jk < zmicedep(ji-1,jj))       ) THEN 
     1231                  mbathy(ji,jj) = 0 ; micedep(ji,jj) = jpk ; icedep(ji,jj) = 0._wp ; bathy(ji,jj) = 0._wp 
     1232               ENDIF 
     1233               ENDIF 
     1234            ENDDO 
     1235        ENDDO 
     1236        WHERE (micedep==jpk) micedep=0  
     1237        IF( lk_mpp ) THEN  
     1238            zbathy(:,:) = FLOAT( micedep(:,:) )  
     1239            CALL lbc_lnk( zbathy, 'T', 1. )  
     1240            micedep(:,:) = INT( zbathy(:,:) )  
     1241            CALL lbc_lnk( icedep, 'T', 1. )  
     1242            CALL lbc_lnk( bathy, 'T', 1. ) 
     1243            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1244            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1245            mbathy(:,:) = INT( zbathy(:,:) ) 
     1246        ENDIF  
     1247        ENDDO 
     1248      END DO  
     1249 
     1250      WHERE (mbathy(:,:) < micedep(:,:)) 
     1251         micedep(:,:) = 0 
     1252         icedep(:,:)  = 0._wp 
     1253         mbathy(:,:)  = 0 
     1254         bathy(:,:)   = 0._wp 
     1255      END WHERE 
     1256 
     1257 
     1258      IF( icompt == 0 ) THEN  
     1259         IF(lwp) WRITE(numout,*)'     no points with ice shelf too close to bathymetry'  
     1260      ELSE  
     1261         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry'  
     1262      ENDIF  
    9081263 
    9091264      ! Scale factors and depth at T- and W-points 
    … …  
    9381293!gm Bug?  check the gdepw_1d 
    9391294                  !       ... on ik 
    940                   gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0   (ji,jj,ik+1) - gdepw_1d(ik) )   & 
    941                      &                           * ((gdept_1d(      ik  ) - gdepw_1d(ik) )   & 
    942                      &                           / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )) 
     1295                  gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   & 
     1296                     &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   & 
     1297                     &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
    9431298                  e3t_0(ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
    9441299                     &                          / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )  
    … …  
    9731328         END DO 
    9741329      END DO 
     1330      ! 
     1331      ! (ISF) Definition of e3t, u, v, w for ISF case 
     1332      DO jj = 1, jpj  
     1333         DO ji = 1, jpi  
     1334            ik = micedep(ji,jj)  
     1335            IF( ik > 1 ) THEN               ! ice shelf point only  
     1336                IF( icedep(ji,jj) < gdepw_1d(ik) )  icedep(ji,jj)= gdepw_1d(ik)  
     1337                gdepw_0(ji,jj,ik) = icedep(ji,jj)  
     1338!gm Bug?  check the gdepw_0  
     1339                !       ... on ik  
     1340                gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   &  
     1341                   &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   &  
     1342                   &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      )  
     1343                e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)  
     1344                e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 
     1345 
     1346                IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only  
     1347                   e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)  
     1348                ENDIF  
     1349                !       ... on ik / ik-1  
     1350                e3w_0  (ji,jj,ik  ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))  
     1351                e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 
     1352! The next line isn't required and doesn't affect results - included for consistency with bathymetry code  
     1353                gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 
     1354             ENDIF  
     1355          END DO  
     1356       END DO  
     1357      !  
     1358      it = 0  
     1359      DO jj = 1, jpj  
     1360         DO ji = 1, jpi  
     1361            ik = micedep(ji,jj)  
     1362            IF( ik > 1 ) THEN               ! ice shelf point only  
     1363               e3tp (ji,jj) = e3t_0(ji,jj,ik  )  
     1364               e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )  
     1365               ! test  
     1366               zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  )  
     1367               IF( zdiff <= 0. .AND. lwp ) THEN   
     1368                  it = it + 1  
     1369                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj  
     1370                  WRITE(numout,*) ' icedep = ', icedep(ji,jj)  
     1371                  WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff  
     1372                  WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj)  
     1373               ENDIF  
     1374            ENDIF  
     1375         END DO  
     1376      END DO  
     1377      ! END (ISF) 
    9751378 
    9761379      ! Scale factors and depth at U-, V-, UW and VW-points 
    … …  
    9911394         END DO 
    9921395      END DO 
     1396      ! (ISF) define e3uw 
     1397      DO jk = 2,jpk                           
     1398         DO jj = 1, jpjm1  
     1399            DO ji = 1, fs_jpim1   ! vector opt.  
     1400! (ISF)** NEEDS CHANGING TO SECOND OPTION FOR ICE SHELF BUT WILL CHANGE RESULTS WITHOUT ICE (DIFFER AT THE 1e-13 LEVEL)  
     1401               e3uw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji+1,jj  ,jk) ) - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji+1,jj  ,jk-1) ) 
     1402               e3vw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji  ,jj+1,jk) ) - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji  ,jj+1,jk-1) ) 
     1403            END DO  
     1404         END DO  
     1405      END DO 
     1406      !End (ISF) 
     1407       
    9931408      CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions 
    9941409      CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 
    … …  
    10321447      
    10331448      ! Compute gdep3w_0 (vertical sum of e3w) 
    1034       gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 
    1035       DO jk = 2, jpk 
    1036          gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk)  
    1037       END DO 
    1038          
     1449      WHERE (micedep == 0) micedep = 1 
     1450      DO jj = 1,jpj 
     1451         DO ji = 1,jpi 
     1452            gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 
     1453            DO jk = 2, micedep(ji,jj) 
     1454               gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
     1455            END DO 
     1456            IF (micedep(ji,jj) .GE. 2) gdep3w_0(ji,jj,micedep(ji,jj)) = icedep(ji,jj) + 0.5_wp * e3w_0(ji,jj,micedep(ji,jj)) 
     1457            DO jk = micedep(ji,jj) + 1, jpk 
     1458               gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
     1459            END DO 
     1460        END DO 
     1461      END DO 
    10391462      !                                               ! ================= ! 
    10401463      IF(lwp .AND. ll_print) THEN                     !   Control print   ! 
    … …  
    10661489      ! 
    10671490      CALL wrk_dealloc( jpi, jpj, jpk, zprt ) 
     1491      CALL wrk_dealloc( jpi, jpj, zmask, zbathy ) 
     1492      CALL wrk_dealloc( jpi, jpj, zmicedep, zmbathy ) 
    10681493      ! 
    10691494      IF( nn_timing == 1 )  CALL timing_stop('zgr_zps') 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90

    r4624 r4666  
    263263                     ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik-1,jp_sal) 
    264264                  ENDIF 
     265                  ik = mikt(ji,jj) 
     266                  IF( ik > 1 ) THEN 
     267                     zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) )  
     268                     ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik+1,jp_tem) 
     269                     ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik+1,jp_sal) 
     270                  END IF 
    265271               END DO 
    266272            END DO 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90

    r4370 r4666  
    7676      ! 
    7777 
    78       IF(lwp) WRITE(numout,*) 
     78      IF(lwp) WRITE(numout,*) ' ' 
    7979      IF(lwp) WRITE(numout,*) 'istate_ini : Initialization of the dynamics and tracers' 
    8080      IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
    … …  
    110110         ELSEIF( cp_cfg == 'gyre' ) THEN          
    111111            CALL istate_gyre                     ! GYRE  configuration : start from pre-defined T-S fields 
     112        ELSEIF( cp_cfg == 'isomip' .OR. cp_cfg == 'isomip2') THEN 
     113            IF(lwp) WRITE(numout,*) 'Initialization of T+S for ISOMIP domain'  
     114            tsn(:,:,:,jp_tem)=-1.9*tmask(:,:,:)          ! ISOMIP configuration : start from constant T+S fields  
     115            tsn(:,:,:,jp_sal)=34.4*tmask(:,:,:) 
     116            tsb(:,:,:,:)=tsn(:,:,:,:)   
    112117         ELSE                                    ! Initial T-S, U-V fields read in files 
    113118            IF ( ln_tsd_init ) THEN              ! read 3D T and S data at nit000 
    … …  
    129134         CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) )        ! before potential and in situ densities 
    130135#if ! defined key_c1d 
    131          IF( ln_zps )   CALL zps_hde( nit000, jpts, tsb, gtsu, gtsv,  & ! zps: before hor. gradient 
    132             &                                       rhd, gru , grv  )   ! of t,s,rd at ocean bottom 
     136         IF( ln_zps )    CALL zps_hde( nit000, jpts, tsb, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient 
     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 
    133139#endif 
    134140         !    

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90

    r3625 r4666  
    5050   REAL(wp), PUBLIC ::   rau0     = 1026._wp         !: volumic mass of reference     [kg/m3] 
    5151#else 
    52    REAL(wp), PUBLIC ::   rau0     = 1035._wp         !: volumic mass of reference     [kg/m3] 
     52   REAL(wp), PUBLIC ::   rau0     = 1028.4_wp       !: volumic mass of reference     [kg/m3] 
    5353#endif 
    5454   REAL(wp), PUBLIC ::   r1_rau0                     !: = 1. / rau0                   [m3/kg] 
Note: See TracChangeset for help on using the changeset viewer.