Changeset 5802


Ignore:
Timestamp:
2015-10-16T19:19:43+02:00 (5 years ago)
Author:
mathiot
Message:

ice sheet coupling: reproducibility fixed in MISOMIP configuration + contrain bathy to be constant during all the run

Location:
branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90

    r5779 r5802  
    200200!      ENDIF 
    201201!!gm end 
    202        
     202      IF ( lk_vvl ) THEN 
     203         IF (lwp) PRINT *, 'cons heat : ', kt, zdiff_hc / zvol_tot, zdiff_sc / zvol_tot  
     204         IF (lwp) PRINT *, 'cons volu : ', kt, zdiff_v2 * 1.e-9       
     205      ELSE 
     206         IF (lwp) PRINT *, 'cons heat : ', kt, zdiff_hc1 * 1.e-20 * rau0 * rcp, zdiff_sc1 * 1.e-9 
     207         IF (lwp) PRINT *, 'cons vol  : ', kt, zdiff_v1 * 1.e-9       
     208      END IF 
    203209      IF( lk_vvl ) THEN 
    204210        CALL iom_put( 'bgtemper' , zdiff_hc / zvol_tot )              ! Temperature variation (C)  
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r5619 r5802  
    253253   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbkt               !: vertical index of the bottom last T- ocean level 
    254254   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbku, mbkv         !: vertical index of the bottom last U- and W- ocean level 
    255    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bathy                     !: ocean depth (meters) 
    256    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_i                   !: interior domain T-point mask 
     255   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bathy              !: ocean depth (meters) 
     256   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_i            !: interior domain T-point mask 
     257   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_h            !: internal domain T-point mask (Figure 8.5 NEMO book) 
    257258   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bmask                     !: land/ocean mask of barotropic stream function 
    258259 
     
    391392 
    392393      ALLOCATE( mbathy(jpi,jpj) , bathy  (jpi,jpj) ,                                       & 
    393          &     tmask_i(jpi,jpj) ,                                                          &  
     394         &     tmask_i(jpi,jpj) , tmask_h(jpi, jpj),                                       &  
    394395         &     ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) , & 
    395396         &     bmask(jpi,jpj)   ,                                                          & 
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90

    r5619 r5802  
    222222      ! -------------------- 
    223223      tmask_i(:,:) = ssmask(:,:)            ! (ISH) tmask_i = 1 even on the ice shelf 
     224 
     225      tmask_h(:,:) = 1._wp                 ! 0 on the halo and 1 elsewhere 
    224226      iif = jpreci                         ! ??? 
    225227      iil = nlci - jpreci + 1 
     
    227229      ijl = nlcj - jprecj + 1 
    228230 
    229       tmask_i( 1 :iif,   :   ) = 0._wp      ! first columns 
    230       tmask_i(iil:jpi,   :   ) = 0._wp      ! last  columns (including mpp extra columns) 
    231       tmask_i(   :   , 1 :ijf) = 0._wp      ! first rows 
    232       tmask_i(   :   ,ijl:jpj) = 0._wp      ! last  rows (including mpp extra rows) 
     231      tmask_h( 1 :iif,   :   ) = 0._wp      ! first columns 
     232      tmask_h(iil:jpi,   :   ) = 0._wp      ! last  columns (including mpp extra columns) 
     233      tmask_h(   :   , 1 :ijf) = 0._wp      ! first rows 
     234      tmask_h(   :   ,ijl:jpj) = 0._wp      ! last  rows (including mpp extra rows) 
    233235 
    234236      ! north fold mask 
     
    241243         IF( mjg(nlej) == jpjglo ) THEN                  ! only half of the nlcj-1 row 
    242244            DO ji = iif+1, iil-1 
    243                tmask_i(ji,nlej-1) = tmask_i(ji,nlej-1) * tpol(mig(ji)) 
     245               tmask_h(ji,nlej-1) = tmask_h(ji,nlej-1) * tpol(mig(ji)) 
    244246            END DO 
    245247         ENDIF 
    246248      ENDIF 
     249      
     250      tmask_i(:,:) = tmask_i(:,:) * tmask_h(:,:) 
     251 
    247252      IF( jperio == 5 .OR. jperio == 6 ) THEN      ! F-point pivot 
    248253         tpol(     1    :jpiglo) = 0._wp 
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r5779 r5802  
    12861286 
    12871287         ! split last cell if possible (only where water column is 2 cell or less) 
    1288          !DO jk = jpkm1, 1, -1 
    1289          !   zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
    1290          !   WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy) 
    1291          !      mbathy(:,:) = jk 
    1292          !      bathy(:,:)  = zmax 
    1293          !   END WHERE 
    1294          !END DO 
     1288         ! if coupled to ice sheet, we do not modify the bathymetry (can be discuss). 
     1289         IF ( .NOT. ln_iscpl) THEN 
     1290            DO jk = jpkm1, 1, -1 
     1291               zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
     1292               WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy) 
     1293                  mbathy(:,:) = jk 
     1294                  bathy(:,:)  = zmax 
     1295               END WHERE 
     1296            END DO 
     1297         END IF 
    12951298  
    12961299         ! split top cell if possible (only where water column is 2 cell or less) 
     
    13101313               ! test bathy 
    13111314               IF (risfdep(ji,jj) .GT. 1) THEN 
    1312                zbathydiff =ABS(bathy(ji,jj)   - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 
    1313                zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  ) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 
     1315                  IF ( .NOT. ln_iscpl ) THEN 
     1316                     zbathydiff  =ABS(bathy(ji,jj)   - (gdepw_1d(mbathy (ji,jj)+1) & 
     1317                         &            + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 
     1318                     zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  ) & 
     1319                         &            - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 
    13141320  
    1315                   IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN 
    1316 !                     IF (zbathydiff .LE. zrisfdepdiff) THEN 
    1317 !                        bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 
    1318 !                        mbathy(ji,jj)= mbathy(ji,jj) + 1 
    1319 !                     ELSE 
     1321                     IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN 
     1322                        IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1323                           bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 
     1324                           mbathy(ji,jj)= mbathy(ji,jj) + 1 
     1325                        ELSE 
     1326                           risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 
     1327                           misfdep(ji,jj) = misfdep(ji,jj) - 1 
     1328                        END IF 
     1329                     ENDIF 
     1330                  ELSE 
     1331                     IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN 
    13201332                        risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 
    13211333                        misfdep(ji,jj) = misfdep(ji,jj) - 1 
    1322 !                     END IF 
     1334                     END IF 
    13231335                  END IF 
    13241336               END IF 
     
    13321344               ! test bathy 
    13331345               IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
    1334 !                  zbathydiff  =ABS(bathy(ji,jj)  - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 
    1335 !                  zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  ) - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 
    1336 !                  IF (zbathydiff .LE. zrisfdepdiff) THEN 
    1337 !                     mbathy(ji,jj) = mbathy(ji,jj) + 1 
    1338 !                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 
    1339 !                  ELSE 
     1346                  IF ( .NOT. ln_iscpl ) THEN  
     1347                     zbathydiff  =ABS(bathy(ji,jj)   - ( gdepw_1d(mbathy (ji,jj)+1) & 
     1348                         &                             + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 
     1349                     zrisfdepdiff=ABS(risfdep(ji,jj) - ( gdepw_1d(misfdep(ji,jj)  ) &  
     1350                         &                             - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 
     1351                     IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1352                        mbathy(ji,jj) = mbathy(ji,jj) + 1 
     1353                        bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 
     1354                     ELSE 
     1355                        misfdep(ji,jj)= misfdep(ji,jj) - 1 
     1356                        risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 
     1357                     END IF 
     1358                  ELSE 
    13401359                     misfdep(ji,jj)= misfdep(ji,jj) - 1 
    1341                      risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 
    1342 !                  END IF 
     1360                     risfdep(ji,jj)= gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 
     1361                  END IF 
    13431362               ENDIF 
    13441363            END DO 
     
    13491368            DO ji = 1, jpim1 
    13501369               IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
    1351 !                  zbathydiff =ABS(bathy(ji,jj    ) - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj  )+1)*e3zps_rat ))) 
    1352 !                  zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 
    1353 !                  IF (zbathydiff .LE. zrisfdepdiff) THEN 
    1354 !                     mbathy(ji,jj) = mbathy(ji,jj) + 1 
    1355 !                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj  )) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj   )+1)*e3zps_rat ) 
    1356 !                  ELSE 
     1370                  IF ( .NOT. ln_iscpl ) THEN  
     1371                     zbathydiff  =ABS(bathy(ji,jj    ) - ( gdepw_1d(mbathy (ji,jj)+1) & 
     1372                          &                              + MIN( e3zps_min, e3t_1d(mbathy (ji,jj  )+1)*e3zps_rat ))) 
     1373                     zrisfdepdiff=ABS(risfdep(ji,jj+1) - ( gdepw_1d(misfdep(ji,jj+1)) & 
     1374                          &                              - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 
     1375                     IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1376                        mbathy(ji,jj) = mbathy(ji,jj) + 1 
     1377                        bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj  )) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj   )+1)*e3zps_rat ) 
     1378                     ELSE 
     1379                        misfdep(ji,jj+1)  = misfdep(ji,jj+1) - 1 
     1380                        risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 
     1381                     END IF 
     1382                  ELSE 
    13571383                     misfdep(ji,jj+1)  = misfdep(ji,jj+1) - 1 
    13581384                     risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 
    1359 !                  END IF 
     1385                  END IF 
    13601386               ENDIF 
    13611387            END DO 
     
    13761402            DO ji = 1, jpim1 
    13771403               IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN 
    1378 !                  zbathydiff =ABS(  bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ))) 
    1379 !                  zrisfdepdiff=ABS(risfdep(ji,jj  ) - (gdepw_1d(misfdep(ji,jj  )  ) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )-1)*e3zps_rat ))) 
    1380 !                  IF (zbathydiff .LE. zrisfdepdiff) THEN 
    1381 !                     mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 
    1382 !                     bathy  (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ) 
    1383 !                  ELSE 
     1404                  IF ( .NOT. ln_iscpl ) THEN  
     1405                     zbathydiff  =ABS(  bathy(ji,jj+1) - ( gdepw_1d(mbathy (ji,jj+1)+1) & 
     1406                           &                             + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ))) 
     1407                     zrisfdepdiff=ABS(risfdep(ji,jj  ) - ( gdepw_1d(misfdep(ji,jj  )  ) & 
     1408                           &                             - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )-1)*e3zps_rat ))) 
     1409                     IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1410                        mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 
     1411                        bathy  (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ) 
     1412                     ELSE 
     1413                        misfdep(ji,jj)   = misfdep(ji,jj) - 1 
     1414                        risfdep(ji,jj)   = gdepw_1d(misfdep(ji,jj  )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )  )*e3zps_rat ) 
     1415                     END IF 
     1416                  ELSE 
    13841417                     misfdep(ji,jj)   = misfdep(ji,jj) - 1 
    13851418                     risfdep(ji,jj)   = gdepw_1d(misfdep(ji,jj  )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )  )*e3zps_rat ) 
    1386 !                  END IF 
     1419                  END IF 
    13871420               ENDIF 
    13881421            END DO 
     
    14051438            DO ji = 1, jpim1 
    14061439               IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
    1407 !                  zbathydiff =ABS(  bathy(ji  ,jj) - (gdepw_1d(mbathy (ji,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji  ,jj)+1)*e3zps_rat ))) 
    1408 !                  zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 
    1409 !                  IF (zbathydiff .LE. zrisfdepdiff) THEN 
    1410 !                     mbathy(ji,jj) = mbathy(ji,jj) + 1 
    1411 !                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 
    1412 !                  ELSE 
     1440                  IF ( .NOT. ln_iscpl ) THEN  
     1441                  zbathydiff  =ABS(  bathy(ji  ,jj) - ( gdepw_1d(mbathy (ji,jj)+1) & 
     1442                       &                              + MIN( e3zps_min, e3t_1d(mbathy (ji  ,jj)+1)*e3zps_rat ))) 
     1443                  zrisfdepdiff=ABS(risfdep(ji+1,jj) - ( gdepw_1d(misfdep(ji+1,jj)) & 
     1444                       &                              - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 
     1445                  IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1446                     mbathy(ji,jj) = mbathy(ji,jj) + 1 
     1447                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 
     1448                  ELSE 
    14131449                     misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 
    14141450                     risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 
    1415 !                  END IF 
     1451                  END IF 
     1452                  ELSE 
     1453                     misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 
     1454                     risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 
     1455                  ENDIF 
    14161456               ENDIF 
    14171457            ENDDO 
     
    14331473            DO ji = 1, jpim1 
    14341474               IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN 
    1435 !                  zbathydiff =ABS(  bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 
    1436 !                  zrisfdepdiff=ABS(risfdep(ji  ,jj) - (gdepw_1d(misfdep(ji  ,jj)  ) - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)-1)*e3zps_rat ))) 
    1437 !                  IF (zbathydiff .LE. zrisfdepdiff) THEN 
    1438 !                     mbathy(ji+1,jj)  = mbathy (ji+1,jj) + 1 
    1439 !                     bathy (ji+1,jj)  = gdepw_1d(mbathy (ji+1,jj)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat ) 
    1440 !                  ELSE 
     1475                  IF ( .NOT. ln_iscpl ) THEN  
     1476                     zbathydiff  =ABS(  bathy(ji+1,jj) - ( gdepw_1d(mbathy (ji+1,jj)+1) & 
     1477                          &                              + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 
     1478                     zrisfdepdiff=ABS(risfdep(ji  ,jj) - ( gdepw_1d(misfdep(ji  ,jj)  ) & 
     1479                          &                              - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)-1)*e3zps_rat ))) 
     1480                     IF (zbathydiff .LE. zrisfdepdiff) THEN 
     1481                        mbathy(ji+1,jj)  = mbathy (ji+1,jj) + 1 
     1482                        bathy (ji+1,jj)  = gdepw_1d(mbathy (ji+1,jj)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat ) 
     1483                     ELSE 
     1484                        misfdep(ji,jj)   = misfdep(ji  ,jj) - 1 
     1485                        risfdep(ji,jj)   = gdepw_1d(misfdep(ji  ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)   )*e3zps_rat ) 
     1486                     END IF 
     1487                  ELSE 
    14411488                     misfdep(ji,jj)   = misfdep(ji  ,jj) - 1 
    14421489                     risfdep(ji,jj)   = gdepw_1d(misfdep(ji  ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)   )*e3zps_rat ) 
    1443 !                  END IF 
     1490                  ENDIF 
    14441491               ENDIF 
    14451492            ENDDO 
     
    15301577               ibtestim1 = zmisfdep(ji-1,jj  ) ; ibtestip1 = zmisfdep(ji+1,jj  ) 
    15311578               ibtestjm1 = zmisfdep(ji  ,jj-1) ; ibtestjp1 = zmisfdep(ji  ,jj+1) 
    1532                IF( zmisfdep(ji,jj) .GE. mbathy(ji-1,jj  ) ) ibtestim1 = jpk!MAX(0, mbathy(ji-1,jj  ) - 1) 
    1533                IF( zmisfdep(ji,jj) .GE. mbathy(ji+1,jj  ) ) ibtestip1 = jpk!MAX(0, mbathy(ji+1,jj  ) - 1) 
    1534                IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj-1) ) ibtestjm1 = jpk!MAX(0, mbathy(ji  ,jj-1) - 1) 
    1535                IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj+1) ) ibtestjp1 = jpk!MAX(0, mbathy(ji  ,jj+1) - 1) 
     1579               IF( zmisfdep(ji,jj) .GE. mbathy(ji-1,jj  ) ) ibtestim1 = jpk 
     1580               IF( zmisfdep(ji,jj) .GE. mbathy(ji+1,jj  ) ) ibtestip1 = jpk 
     1581               IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj-1) ) ibtestjm1 = jpk 
     1582               IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj+1) ) ibtestjp1 = jpk 
    15361583               ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 
    15371584               IF( ibtest == jpk .AND. misfdep(ji,jj) .GE. 2) THEN 
     
    15621609               ibtestim1 = zmbathy(ji-1,jj  ) ; ibtestip1 = zmbathy(ji+1,jj  ) 
    15631610               ibtestjm1 = zmbathy(ji  ,jj-1) ; ibtestjp1 = zmbathy(ji  ,jj+1) 
    1564                IF( zmbathy(ji,jj) .LT. misfdep(ji-1,jj  ) ) ibtestim1 = 0!MIN(jpk-1, misfdep(ji-1,jj  ) + 1) 
     1611               IF( zmbathy(ji,jj) .LT. misfdep(ji-1,jj  ) ) ibtestim1 = 0 
    15651612               IF( zmbathy(ji,jj) .LT. misfdep(ji+1,jj  ) ) ibtestip1 = 0 
    15661613               IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj-1) ) ibtestjm1 = 0 
     
    17561803               ENDIF  
    17571804            !       ... on ik / ik-1  
    1758                e3w_0  (ji,jj,ik  ) = e3t_0  (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))  
     1805               e3w_0  (ji,jj,ik  ) = e3t_0  (ji,jj,ik)  
    17591806               e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 
    17601807! The next line isn't required and doesn't affect results - included for consistency with bathymetry code  
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90

    r5790 r5802  
    6464      REAL(wp):: zjip1_ratio, zjim1_ratio, zjjp1_ratio, zjjm1_ratio 
    6565      !! 
    66       REAL(wp), DIMENSION(:,:    ), POINTER :: zde3t 
     66      REAL(wp), DIMENSION(:,:    ), POINTER :: zde3t, zdtem, zdsal 
    6767      REAL(wp), DIMENSION(:,:    ), POINTER :: zssh0   
    6868      REAL(wp), DIMENSION(:,:,:  ), POINTER :: ztmp3d 
     
    7474 
    7575      CALL wrk_alloc(jpi,jpj,jpk,   ztmp3d ) 
    76       CALL wrk_alloc(jpi,jpj,       zde3t ) 
     76      CALL wrk_alloc(jpi,jpj,       zde3t , zdtem, zdsal ) 
    7777      CALL wrk_alloc(jpi,jpj,       zssh0  ) 
    7878 
    7979    ! get unbalance (volume heat and salt) 
    8080    ! initialisation 
    81        zde3t   (:,:)     = 0.0_wp 
    82        pvol_flx(:,:,:  ) = 0.0_wp 
    83        pts_flx (:,:,:,:) = 0.0_wp 
    84  
    85        zsum = glob_sum_full(pvol_flx(:,:,:)      ) * rn_fiscpl * rn_rdt 
    86        IF (lwp) PRINT *, 'total volume correction 0 = ',zsum 
    87        zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt 
    88        IF (lwp) PRINT *, 'total heat correction 0 = ',zsum 
    89        zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt 
    90        IF (lwp) PRINT *, 'total salt correction 0 = ',zsum 
    91  
    92        ! mask tsn and tsb (should be useless) 
    93        tsb(:,:,:,jp_tem)=tsb(:,:,:,jp_tem)*ptmask_b(:,:,:); tsn(:,:,:,jp_tem)=tsn(:,:,:,jp_tem)*tmask(:,:,:); 
    94        tsb(:,:,:,jp_sal)=tsb(:,:,:,jp_sal)*ptmask_b(:,:,:); tsn(:,:,:,jp_sal)=tsn(:,:,:,jp_sal)*tmask(:,:,:); 
    95  
    96        ! diagnose non conservation of heat, salt and volume  
    97        r1_tiscpl = 1._wp / (prdt_iscpl * rn_rdt)  
    98        zssh0(:,:)        = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:) 
    99        IF ( lk_vvl ) zssh0 = 0.0_wp 
    100        DO jk = 1,jpk-1 
    101           DO ji = 2,jpi-1 
    102              DO jj = 2,jpj-1 
    103                 ! volume differences 
    104                 zde3t(ji,jj) = fse3t_n(ji,jj,jk) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk); 
    105                  
    106                 ! shh changes 
    107                 IF ( ptmask_b(ji,jj,jk) == 1 .OR. tmask(ji,jj,jk) == 1 ) THEN  
    108                    zde3t(ji,jj) = zde3t(ji,jj) + zssh0(ji,jj) 
    109                    zssh0(ji,jj) = 0._wp 
    110                 END IF 
    111  
    112                 ! ocean cell now 
    113                 ! case where we open, enlarge or thin a cell : 
    114                 pvol_flx(ji,jj,jk)       =                          zde3t(ji,jj) * r1_tiscpl 
    115                 pts_flx (ji,jj,jk,jp_sal)=   tsn(ji,jj,jk,jp_sal) * zde3t(ji,jj) * r1_tiscpl  
    116                 pts_flx (ji,jj,jk,jp_tem)=   tsn(ji,jj,jk,jp_tem) * zde3t(ji,jj) * r1_tiscpl 
    117              END DO 
    118           END DO 
    119        END DO 
    120        ! glob_sum_full because with glob summ some data can be masked. WARNING the halo have to be set at 0 
    121        PRINT *, 'test ', narea, SUM(pvol_flx(:,:,:)) * rn_fiscpl * rn_rdt, SUM(pvol_flx(2:jpi-1,2:jpj-1,:)) * rn_fiscpl * rn_rdt 
    122        zsum = glob_sum_full(pvol_flx(:,:,:)      ) * rn_fiscpl * rn_rdt 
    123        IF (lwp) PRINT *, 'total volume correction 1 = ',zsum 
    124        zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt 
    125        IF (lwp) PRINT *, 'total heat correction 1 = ',zsum 
    126        zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt 
    127        IF (lwp) PRINT *, 'total salt correction 1 = ',zsum 
    128  
    129        zssh0(:,:)        = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:) 
    130        IF ( lk_vvl ) zssh0 = 0.0_wp 
    131        DO jk = 1,jpk-1 
    132           DO ji = 2,jpi-1 
    133              DO jj = 2,jpj-1 
    134                 ! volume differences 
    135                 zde3t(ji,jj) = fse3t_n(ji,jj,jk) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk); 
    136                  
    137                 ! shh changes 
    138                 IF ( ptmask_b(ji,jj,jk) == 1 .OR. tmask(ji,jj,jk) == 1 ) THEN  
    139                    zde3t(ji,jj) = zde3t(ji,jj) + zssh0(ji,jj) 
    140                    zssh0(ji,jj) = 0._wp 
    141                 END IF 
    142  
    143                 ! ocean cell before and mask cell now 
    144                 IF ( tmask(ji,jj,jk) == 0._wp .AND. ptmask_b(ji,jj,jk) == 1._wp ) THEN 
    145                    ! case where we close a cell and adjacent cell open 
    146                    pvol_flx(ji,jj,jk)       =                         zde3t(ji,jj) * r1_tiscpl 
    147                    pts_flx (ji,jj,jk,jp_sal)=  tsb(ji,jj,jk,jp_sal) * zde3t(ji,jj) * r1_tiscpl  
    148                    pts_flx (ji,jj,jk,jp_tem)=  tsb(ji,jj,jk,jp_tem) * zde3t(ji,jj) * r1_tiscpl 
    149  
    150                    jip1=ji+1 ; jim1=ji-1 ; jjp1=jj+1 ; jjm1=jj-1 ; 
    151  
    152                    zsum =   e12t(ji  ,jjp1) * tmask(ji  ,jjp1,jk) + e12t(ji  ,jjm1) * tmask(ji  ,jjm1,jk) & 
    153                      &    + e12t(jim1,jj  ) * tmask(jim1,jj  ,jk) + e12t(jip1,jj  ) * tmask(jip1,jj  ,jk)  
    154  
    155                    IF ( zsum .NE. 0._wp ) THEN 
    156                       zjip1_ratio = e12t(jip1,jj  ) * tmask(jip1,jj  ,jk) / zsum 
    157                       zjim1_ratio = e12t(jim1,jj  ) * tmask(jim1,jj  ,jk) / zsum 
    158                       zjjp1_ratio = e12t(ji  ,jjp1) * tmask(ji  ,jjp1,jk) / zsum 
    159                       zjjm1_ratio = e12t(ji  ,jjm1) * tmask(ji  ,jjm1,jk) / zsum 
    160  
    161                       pvol_flx(ji  ,jjp1,jk       ) = pvol_flx(ji  ,jjp1,jk       ) + pvol_flx(ji,jj,jk       ) * zjjp1_ratio 
    162                       pvol_flx(ji  ,jjm1,jk       ) = pvol_flx(ji  ,jjm1,jk       ) + pvol_flx(ji,jj,jk       ) * zjjm1_ratio 
    163                       pvol_flx(jip1,jj  ,jk       ) = pvol_flx(jip1,jj  ,jk       ) + pvol_flx(ji,jj,jk       ) * zjip1_ratio 
    164                       pvol_flx(jim1,jj  ,jk       ) = pvol_flx(jim1,jj  ,jk       ) + pvol_flx(ji,jj,jk       ) * zjim1_ratio 
    165                       pts_flx (ji  ,jjp1,jk,jp_sal) = pts_flx (ji  ,jjp1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjp1_ratio 
    166                       pts_flx (ji  ,jjm1,jk,jp_sal) = pts_flx (ji  ,jjm1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjm1_ratio 
    167                       pts_flx (jip1,jj  ,jk,jp_sal) = pts_flx (jip1,jj  ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjip1_ratio 
    168                       pts_flx (jim1,jj  ,jk,jp_sal) = pts_flx (jim1,jj  ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjim1_ratio 
    169                       pts_flx (ji  ,jjp1,jk,jp_tem) = pts_flx (ji  ,jjp1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjp1_ratio 
    170                       pts_flx (ji  ,jjm1,jk,jp_tem) = pts_flx (ji  ,jjm1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjm1_ratio 
    171                       pts_flx (jip1,jj  ,jk,jp_tem) = pts_flx (jip1,jj  ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjip1_ratio 
    172                       pts_flx (jim1,jj  ,jk,jp_tem) = pts_flx (jim1,jj  ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjim1_ratio 
    173  
    174                       ! set to 0 the cell we distributed over neigbourg cells 
    175                       pvol_flx(ji,jj,jk       ) = 0._wp 
    176                       pts_flx (ji,jj,jk,jp_sal) = 0._wp 
    177                       pts_flx (ji,jj,jk,jp_tem) = 0._wp 
    178  
    179                    ELSE IF (zsum .EQ. 0._wp ) THEN 
    180                       ! case where we close a cell and no adjacent cell open 
    181                       ! check if the cell beneath is wet 
    182                       IF ( tmask(ji,jj,jk+1) .EQ. 1._wp ) THEN 
    183                          pvol_flx(ji,jj,jk+1)       =  pvol_flx(ji,jj,jk+1)        + pvol_flx(ji,jj,jk) 
    184                          pts_flx (ji,jj,jk+1,jp_sal)=  pts_flx (ji,jj,jk+1,jp_sal) + pts_flx (ji,jj,jk,jp_sal) 
    185                          pts_flx (ji,jj,jk+1,jp_tem)=  pts_flx (ji,jj,jk+1,jp_tem) + pts_flx (ji,jj,jk,jp_tem) 
    186  
    187                          ! set to 0 the cell we distributed over neigbourg cells 
    188                          pvol_flx(ji,jj,jk       ) = 0._wp 
    189                          pts_flx (ji,jj,jk,jp_sal) = 0._wp 
    190                          pts_flx (ji,jj,jk,jp_tem) = 0._wp 
    191                       ELSE 
    192                       ! case no adjacent cell on the horizontal and on the vertical 
    193                          PRINT *,        'W A R N I N G iscpl: no adjacent cell on the vertical and horizontal' 
    194                          PRINT *,        '                     ',mig(ji),' ',mjg(jj),' ',jk 
    195                          PRINT *,        '                     ',ji,' ',jj,' ',jk,' ',narea 
    196                          PRINT *,        ' we are now looking for the closest wet cell on the horizontal ' 
    197                       ! We deal with this points later. 
    198                       END IF 
    199                    END IF 
    200                 END IF 
    201              END DO 
    202           END DO 
    203        END DO 
    204  
    205        zsum = glob_sum_full(pvol_flx(:,:,:)      ) * rn_fiscpl * rn_rdt 
    206        IF (lwp) PRINT *, 'total volume correction 2 = ',zsum 
    207        zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt 
    208        IF (lwp) PRINT *, 'total heat correction 2 = ',zsum 
    209        zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt 
    210        IF (lwp) PRINT *, 'total salt correction 2 = ',zsum 
    211  
    212        ! allocation and initialisation of the list of problematic point 
    213        ALLOCATE(vnpts(jpnij)) 
    214        vnpts(:)=0 
    215  
    216        ! fill narea location with the number of problematic point 
    217        DO jk = 1,jpk-1 
    218           DO ji = 2,jpi-1 
    219              DO jj = 2,jpj-1 
    220                 IF (          ptmask_b(ji,jj,jk)      == 1 .AND. SUM(tmask(ji-1:ji+1,jj,jk)) == 0 & 
    221                 &   .AND. SUM(tmask(ji,jj-1:jj+1,jk)) == 0 .AND.     tmask(ji,jj,jk+1)       == 0 ) THEN 
    222                    vnpts(narea) = vnpts(narea) + 1  
    223                 END IF 
    224              END DO 
    225           END DO 
    226        END DO 
    227  
    228        ! build array of total problematic point on each cpu (share to each cpu) 
    229        CALL mpp_max(vnpts,jpnij)  
    230  
    231        ! size of the new variable 
    232        npts  = SUM(vnpts)     
    233         
    234        ! allocation of the coordinates, correction, index vector for the problematic points  
    235        ALLOCATE(ixpts(npts), iypts(npts), izpts(npts), zcorr_vol(npts), zcorr_sal(npts), zcorr_tem(npts), zlon(npts), zlat(npts)) 
    236        ixpts(:) = -9999 ; iypts(:) = -9999 ; izpts(:) = -9999 ; zlon(:) = -1.0e20 ; zlat(:) = -1.0e20 
    237        zcorr_vol(:) = 0.0_wp 
    238        zcorr_sal(:) = 0.0_wp 
    239        zcorr_tem(:) = 0.0_wp 
    240  
    241        ! fill new variable 
    242        jpts = SUM(vnpts(1:narea-1)) 
    243        DO jk = 1,jpk-1 
    244           DO ji = 2,jpi-1 
    245              DO jj = 2,jpj-1 
    246                 IF (          ptmask_b(ji,jj,jk)      == 1 .AND. SUM(tmask(ji-1:ji+1,jj,jk)) == 0 & 
    247                 &   .AND. SUM(tmask(ji,jj-1:jj+1,jk)) == 0 .AND.     tmask(ji,jj,jk+1)       == 0 ) THEN 
    248                    jpts = jpts + 1  ! positioning in the vnpts vector for the area narea 
    249                    PRINT *, 'corrected point ', narea, ji, jj, jk, jpts 
    250                    ixpts(jpts) = ji           ; iypts(jpts) = jj ; izpts(jpts) = jk 
    251                    zlon (jpts) = glamt(ji,jj) ; zlat (jpts) = gphit(ji,jj) 
    252                    zcorr_vol(jpts) = pvol_flx(ji,jj,jk) 
    253                    zcorr_sal(jpts) = pts_flx (ji,jj,jk,jp_sal) 
    254                    zcorr_tem(jpts) = pts_flx (ji,jj,jk,jp_tem) 
    255                    ! set flx to 0 (safer) 
    256                    pvol_flx(ji,jj,jk       ) = 0.0_wp 
    257                    pts_flx (ji,jj,jk,jp_sal) = 0.0_wp 
    258                    pts_flx (ji,jj,jk,jp_tem) = 0.0_wp 
    259                    PRINT *, zcorr_vol(jpts)*rn_fiscpl*rn_rdt, zcorr_sal(jpts)*rn_fiscpl*rn_rdt, zcorr_tem(jpts)*rn_fiscpl*rn_rdt 
    260                 END IF 
    261              END DO 
    262           END DO 
    263        END DO 
    264  
    265        ! build array of total problematic point on each cpu (share to each cpu) 
    266        CALL mpp_max(zlat ,npts) 
    267        CALL mpp_max(zlon ,npts) 
    268        CALL mpp_max(izpts,npts) 
    269  
    270        ! put correction term in the closest cell           
    271        PRINT *, 'corrected point1 ', narea, zlon, zlat, izpts 
    272        DO jpts = 1,npts 
    273           CALL dom_ngb(zlon(jpts), zlat(jpts), ixpts(jpts), iypts(jpts),'T', izpts(jpts)) 
    274           PRINT *, 'corrected point2 ', narea, jpts, ixpts(jpts), iypts(jpts), izpts(jpts) 
    275           DO jj = mj0(iypts(jpts)),mj1(iypts(jpts)) 
    276              DO ji = mi0(ixpts(jpts)),mi1(ixpts(jpts)) 
    277                 jk = izpts(jpts) 
    278                 pvol_flx(ji,jj,jk)        =  pvol_flx(ji,jj,jk       ) + zcorr_vol(jpts) 
    279                 pts_flx (ji,jj,jk,jp_sal) =  pts_flx (ji,jj,jk,jp_sal) + zcorr_sal(jpts) 
    280                 pts_flx (ji,jj,jk,jp_tem) =  pts_flx (ji,jj,jk,jp_tem) + zcorr_tem(jpts) 
    281              END DO 
    282           END DO 
    283        END DO 
    284        ! deallocate variables  
    285        DEALLOCATE(vnpts) 
    286        DEALLOCATE(ixpts, iypts, izpts, zcorr_vol, zcorr_sal, zcorr_tem, zlon, zlat) 
    287       
    288        ! add contribution store on the hallo (lbclnk remove one of the contribution) 
    289        pvol_flx(:,:,:       ) = pvol_flx(:,:,:       ) * tmask(:,:,:) 
    290        pts_flx (:,:,:,jp_tem) = pts_flx (:,:,:,jp_tem) * tmask(:,:,:) 
    291        pts_flx (:,:,:,jp_sal) = pts_flx (:,:,:,jp_sal) * tmask(:,:,:) 
    292  
    293        CALL lbc_sum(pvol_flx(:,:,:       ),'T',1.) 
    294        CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1.) 
    295        CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1.) 
    296  
    297        ! CHECK vol !!!!!!!!! warning tmask_i wrong if deals with before value, so glob_sum wrong for before value!!!! 
    298        zsumn = glob_sum     ( fse3t_n(:,:,:) * tmask  (:,:,:)) - glob_sum(pvol_flx(:,:,:)) * rn_fiscpl * rn_rdt 
    299        ztmp3d(:,:,:) = 0.0 
    300        ztmp3d(2:jpi-1,2:jpj-1,:) = pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:) 
    301        zsumb = glob_sum_full(ztmp3d)  
    302        zsum  = glob_sum     ( pvol_flx(:,:,:) * rn_fiscpl * rn_rdt) 
    303        IF (lwp) PRINT *, 'CHECK vol = ',zsumn, zsumb, zsumn - zsumb, zsum 
    304        ! CHECK salt 
    305        zsumn = glob_sum( tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) *  tmask  (:,:,:)) - glob_sum(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt 
    306        zsumb = glob_sum( tsb(:,:,:,jp_sal) *  pe3t_b(:,:,:) * ptmask_b(:,:,:)) 
    307        zsum  = glob_sum( pts_flx(:,:,:,jp_sal)*rn_fiscpl * rn_rdt) 
    308        IF (lwp) PRINT *, 'CHECK salt = ',zsumn, zsumb, zsumn - zsumb, zsum 
    309        ! CHECK heat 
    310        zsumn = glob_sum( tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) *  tmask  (:,:,:)) - glob_sum(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt 
    311        zsumb = glob_sum( tsb(:,:,:,jp_tem) *  pe3t_b(:,:,:) * ptmask_b(:,:,:)) 
    312        zsum  = glob_sum( pts_flx(:,:,:,jp_tem)*rn_fiscpl * rn_rdt) 
    313        IF (lwp) PRINT *, 'CHECK heat = ',zsumn, zsumb, zsumn - zsumb, zsum 
    314     !!  
    315     CALL wrk_dealloc(jpi,jpj,jpk,   ztmp3d )  
    316     CALL wrk_dealloc(jpi,jpj,       zde3t  )  
    317     CALL wrk_dealloc(jpi,jpj,       zssh0  )  
     81      zde3t   (:,:)     = 0.0_wp 
     82      pvol_flx(:,:,:  ) = 0.0_wp 
     83      pts_flx (:,:,:,:) = 0.0_wp 
     84 
     85      ! mask tsn and tsb (should be useless) 
     86      tsb(:,:,:,jp_tem)=tsb(:,:,:,jp_tem)*ptmask_b(:,:,:); tsn(:,:,:,jp_tem)=tsn(:,:,:,jp_tem)*tmask(:,:,:); 
     87      tsb(:,:,:,jp_sal)=tsb(:,:,:,jp_sal)*ptmask_b(:,:,:); tsn(:,:,:,jp_sal)=tsn(:,:,:,jp_sal)*tmask(:,:,:); 
     88 
     89      ! diagnose non conservation of heat, salt and volume  
     90      r1_tiscpl = 1._wp / (prdt_iscpl * rn_rdt)  
     91 
     92      zssh0(:,:)        = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:) 
     93      IF ( lk_vvl ) zssh0 = 0.0_wp ! already include in the levels by definition 
     94 
     95      DO jk = 1,jpk-1 
     96         DO ji = 2,jpi-1 
     97            DO jj = 2,jpj-1 
     98               IF (tmask_h(ji,jj) == 1._wp) THEN 
     99 
     100               ! volume differences 
     101               zde3t(ji,jj) = fse3t_n(ji,jj,jk) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk) 
     102 
     103               ! heat diff 
     104               zdtem(ji,jj) = tsn(ji,jj,jk,jp_tem) * fse3t_n(ji,jj,jk) * tmask(ji,jj,jk)   & 
     105                            - tsb(ji,jj,jk,jp_tem) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk) 
     106               ! salt diff 
     107               zdsal(ji,jj) = tsn(ji,jj,jk,jp_sal) * fse3t_n(ji,jj,jk) * tmask(ji,jj,jk)   & 
     108                            - tsb(ji,jj,jk,jp_sal) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk) 
     109                
     110               ! shh changes 
     111               IF ( ptmask_b(ji,jj,jk) == 1 .OR. tmask(ji,jj,jk) == 1 ) THEN  
     112                  zde3t(ji,jj) = zde3t(ji,jj) + zssh0(ji,jj) 
     113                  zssh0(ji,jj) = 0._wp 
     114               END IF 
     115 
     116               ! volume, heat and salt differences in each cell  
     117               pvol_flx(ji,jj,jk)       =   pvol_flx(ji,jj,jk)        + zde3t(ji,jj) * r1_tiscpl 
     118               pts_flx (ji,jj,jk,jp_sal)=   pts_flx (ji,jj,jk,jp_sal) + zdsal(ji,jj) * r1_tiscpl  
     119               pts_flx (ji,jj,jk,jp_tem)=   pts_flx (ji,jj,jk,jp_tem) + zdtem(ji,jj) * r1_tiscpl 
     120 
     121               IF ( tmask(ji,jj,jk) == 0._wp .AND. ptmask_b(ji,jj,jk) == 1._wp ) THEN 
     122                  ! case where we close a cell: check if the neighbour cells are wet  
     123 
     124                  jip1=ji+1 ; jim1=ji-1 ; jjp1=jj+1 ; jjm1=jj-1 ; 
     125 
     126                  zsum =   e12t(ji  ,jjp1) * tmask(ji  ,jjp1,jk) + e12t(ji  ,jjm1) * tmask(ji  ,jjm1,jk) & 
     127                    &    + e12t(jim1,jj  ) * tmask(jim1,jj  ,jk) + e12t(jip1,jj  ) * tmask(jip1,jj  ,jk)  
     128 
     129                  IF ( zsum .NE. 0._wp ) THEN 
     130                     zjip1_ratio = e12t(jip1,jj  ) * tmask(jip1,jj  ,jk) / zsum 
     131                     zjim1_ratio = e12t(jim1,jj  ) * tmask(jim1,jj  ,jk) / zsum 
     132                     zjjp1_ratio = e12t(ji  ,jjp1) * tmask(ji  ,jjp1,jk) / zsum 
     133                     zjjm1_ratio = e12t(ji  ,jjm1) * tmask(ji  ,jjm1,jk) / zsum 
     134 
     135                     pvol_flx(ji  ,jjp1,jk       ) = pvol_flx(ji  ,jjp1,jk       ) + pvol_flx(ji,jj,jk       ) * zjjp1_ratio 
     136                     pvol_flx(ji  ,jjm1,jk       ) = pvol_flx(ji  ,jjm1,jk       ) + pvol_flx(ji,jj,jk       ) * zjjm1_ratio 
     137                     pvol_flx(jip1,jj  ,jk       ) = pvol_flx(jip1,jj  ,jk       ) + pvol_flx(ji,jj,jk       ) * zjip1_ratio 
     138                     pvol_flx(jim1,jj  ,jk       ) = pvol_flx(jim1,jj  ,jk       ) + pvol_flx(ji,jj,jk       ) * zjim1_ratio 
     139                     pts_flx (ji  ,jjp1,jk,jp_sal) = pts_flx (ji  ,jjp1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjp1_ratio 
     140                     pts_flx (ji  ,jjm1,jk,jp_sal) = pts_flx (ji  ,jjm1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjm1_ratio 
     141                     pts_flx (jip1,jj  ,jk,jp_sal) = pts_flx (jip1,jj  ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjip1_ratio 
     142                     pts_flx (jim1,jj  ,jk,jp_sal) = pts_flx (jim1,jj  ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjim1_ratio 
     143                     pts_flx (ji  ,jjp1,jk,jp_tem) = pts_flx (ji  ,jjp1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjp1_ratio 
     144                     pts_flx (ji  ,jjm1,jk,jp_tem) = pts_flx (ji  ,jjm1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjm1_ratio 
     145                     pts_flx (jip1,jj  ,jk,jp_tem) = pts_flx (jip1,jj  ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjip1_ratio 
     146                     pts_flx (jim1,jj  ,jk,jp_tem) = pts_flx (jim1,jj  ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjim1_ratio 
     147 
     148                     ! set to 0 the cell we distributed over neigbourg cells 
     149                     pvol_flx(ji,jj,jk       ) = 0._wp 
     150                     pts_flx (ji,jj,jk,jp_sal) = 0._wp 
     151                     pts_flx (ji,jj,jk,jp_tem) = 0._wp 
     152 
     153                  ELSE IF (zsum .EQ. 0._wp ) THEN 
     154                     ! case where we close a cell and no adjacent cell open 
     155                     ! check if the cell beneath is wet 
     156                     IF ( tmask(ji,jj,jk+1) .EQ. 1._wp ) THEN 
     157                        pvol_flx(ji,jj,jk+1)       =  pvol_flx(ji,jj,jk+1)        + pvol_flx(ji,jj,jk) 
     158                        pts_flx (ji,jj,jk+1,jp_sal)=  pts_flx (ji,jj,jk+1,jp_sal) + pts_flx (ji,jj,jk,jp_sal) 
     159                        pts_flx (ji,jj,jk+1,jp_tem)=  pts_flx (ji,jj,jk+1,jp_tem) + pts_flx (ji,jj,jk,jp_tem) 
     160 
     161                        ! set to 0 the cell we distributed over neigbourg cells 
     162                        pvol_flx(ji,jj,jk       ) = 0._wp 
     163                        pts_flx (ji,jj,jk,jp_sal) = 0._wp 
     164                        pts_flx (ji,jj,jk,jp_tem) = 0._wp 
     165                     ELSE 
     166                     ! case no adjacent cell on the horizontal and on the vertical 
     167                        PRINT *,        'W A R N I N G iscpl: no adjacent cell on the vertical and horizontal' 
     168                        PRINT *,        '                     ',mig(ji),' ',mjg(jj),' ',jk 
     169                        PRINT *,        '                     ',ji,' ',jj,' ',jk,' ',narea 
     170                        PRINT *,        ' we are now looking for the closest wet cell on the horizontal ' 
     171                     ! We deal with this points later. 
     172                     END IF 
     173                  END IF 
     174               END IF 
     175               END IF 
     176            END DO 
     177         END DO 
     178      END DO 
     179 
     180      CALL lbc_sum(pvol_flx(:,:,:       ),'T',1.) 
     181      CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1.) 
     182      CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1.) 
     183 
     184      zsum = glob_sum_full(pvol_flx(:,:,:)      ) * rn_fiscpl * rn_rdt 
     185      IF (lwp) PRINT *, 'total volume correction 21 = ',zsum 
     186      zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt 
     187      IF (lwp) PRINT *, 'total heat correction 21 = ',zsum 
     188      zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt 
     189      IF (lwp) PRINT *, 'total salt correction 21 = ',zsum 
     190 
     191      ! if no neighbour wet cell in case of 2close a cell", need to find the nearest wet point  
     192      ! allocation and initialisation of the list of problematic point 
     193      ALLOCATE(vnpts(jpnij)) 
     194      vnpts(:)=0 
     195 
     196      ! fill narea location with the number of problematic point 
     197      DO jk = 1,jpk-1 
     198         DO ji = 2,jpi-1 
     199            DO jj = 2,jpj-1 
     200               IF (          ptmask_b(ji,jj,jk)      == 1 .AND. SUM(tmask(ji-1:ji+1,jj,jk)) == 0 & 
     201               &   .AND. SUM(tmask(ji,jj-1:jj+1,jk)) == 0 .AND.     tmask(ji,jj,jk+1)       == 0 & 
     202               &   .AND. tmask_h(ji,jj) == 1._wp ) THEN 
     203                  vnpts(narea) = vnpts(narea) + 1  
     204               END IF 
     205            END DO 
     206         END DO 
     207      END DO 
     208 
     209      ! build array of total problematic point on each cpu (share to each cpu) 
     210      CALL mpp_max(vnpts,jpnij)  
     211 
     212      ! size of the new variable 
     213      npts  = SUM(vnpts)     
     214       
     215      ! allocation of the coordinates, correction, index vector for the problematic points  
     216      ALLOCATE(ixpts(npts), iypts(npts), izpts(npts), zcorr_vol(npts), zcorr_sal(npts), zcorr_tem(npts), zlon(npts), zlat(npts)) 
     217      ixpts(:) = -9999 ; iypts(:) = -9999 ; izpts(:) = -9999 ; zlon(:) = -1.0e20 ; zlat(:) = -1.0e20 
     218      zcorr_vol(:) = 0.0_wp 
     219      zcorr_sal(:) = 0.0_wp 
     220      zcorr_tem(:) = 0.0_wp 
     221 
     222      ! fill new variable 
     223      jpts = SUM(vnpts(1:narea-1)) 
     224      DO jk = 1,jpk-1 
     225         DO ji = 2,jpi-1 
     226            DO jj = 2,jpj-1 
     227               IF (          ptmask_b(ji,jj,jk)      == 1 .AND. SUM(tmask(ji-1:ji+1,jj,jk)) == 0 & 
     228               &   .AND. SUM(tmask(ji,jj-1:jj+1,jk)) == 0 .AND.     tmask(ji,jj,jk+1)       == 0 & 
     229               &   .AND. tmask_h(ji,jj) == 1 ) THEN 
     230                  jpts = jpts + 1  ! positioning in the vnpts vector for the area narea 
     231                  PRINT *, 'corrected point ', narea, ji, jj, jk, jpts 
     232                  ixpts(jpts) = ji           ; iypts(jpts) = jj ; izpts(jpts) = jk 
     233                  zlon (jpts) = glamt(ji,jj) ; zlat (jpts) = gphit(ji,jj) 
     234                  zcorr_vol(jpts) = pvol_flx(ji,jj,jk) 
     235                  zcorr_sal(jpts) = pts_flx (ji,jj,jk,jp_sal) 
     236                  zcorr_tem(jpts) = pts_flx (ji,jj,jk,jp_tem) 
     237                  ! set flx to 0 (safer) 
     238                  pvol_flx(ji,jj,jk       ) = 0.0_wp 
     239                  pts_flx (ji,jj,jk,jp_sal) = 0.0_wp 
     240                  pts_flx (ji,jj,jk,jp_tem) = 0.0_wp 
     241                  PRINT *, zcorr_vol(jpts)*rn_fiscpl*rn_rdt, zcorr_sal(jpts)*rn_fiscpl*rn_rdt, zcorr_tem(jpts)*rn_fiscpl*rn_rdt 
     242               END IF 
     243            END DO 
     244         END DO 
     245      END DO 
     246 
     247      ! build array of total problematic point on each cpu (share to each cpu) 
     248      CALL mpp_max(zlat ,npts) 
     249      CALL mpp_max(zlon ,npts) 
     250      CALL mpp_max(izpts,npts) 
     251 
     252      ! put correction term in the closest cell           
     253      PRINT *, 'corrected point1 ', narea, zlon, zlat, izpts 
     254      DO jpts = 1,npts 
     255         CALL dom_ngb(zlon(jpts), zlat(jpts), ixpts(jpts), iypts(jpts),'T', izpts(jpts)) 
     256         PRINT *, 'corrected point2 ', narea, jpts, ixpts(jpts), iypts(jpts), izpts(jpts) 
     257         DO jj = mj0(iypts(jpts)),mj1(iypts(jpts)) 
     258            DO ji = mi0(ixpts(jpts)),mi1(ixpts(jpts)) 
     259               jk = izpts(jpts) 
     260               pvol_flx(ji,jj,jk)        =  pvol_flx(ji,jj,jk       ) + zcorr_vol(jpts) 
     261               pts_flx (ji,jj,jk,jp_sal) =  pts_flx (ji,jj,jk,jp_sal) + zcorr_sal(jpts) 
     262               pts_flx (ji,jj,jk,jp_tem) =  pts_flx (ji,jj,jk,jp_tem) + zcorr_tem(jpts) 
     263            END DO 
     264         END DO 
     265      END DO 
     266      ! deallocate variables  
     267      DEALLOCATE(vnpts) 
     268      DEALLOCATE(ixpts, iypts, izpts, zcorr_vol, zcorr_sal, zcorr_tem, zlon, zlat) 
     269     
     270      ! add contribution store on the hallo (lbclnk remove one of the contribution) 
     271      pvol_flx(:,:,:       ) = pvol_flx(:,:,:       ) * tmask(:,:,:) 
     272      pts_flx (:,:,:,jp_sal) = pts_flx (:,:,:,jp_sal) * tmask(:,:,:) 
     273      pts_flx (:,:,:,jp_tem) = pts_flx (:,:,:,jp_tem) * tmask(:,:,:) 
     274 
     275      CALL lbc_sum(pvol_flx(:,:,:       ),'T',1.) 
     276      CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1.) 
     277      CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1.) 
     278 
     279      ! CHECK vol !!!!!!!!! warning tmask_i wrong if deals with before value, so glob_sum wrong for before value!!!! 
     280      zsum  = glob_sum_full( pvol_flx(:,:,:) ) 
     281      IF (lwp) PRINT *, 'CHECK vol = ',zsum 
     282      ! CHECK salt 
     283      zsum  = glob_sum( pts_flx(:,:,:,jp_sal) )  
     284      IF (lwp) PRINT *, 'CHECK salt = ',zsum 
     285      ! CHECK heat 
     286      zsum  = glob_sum( pts_flx(:,:,:,jp_tem) ) 
     287      IF (lwp) PRINT *, 'CHECK heat = ',zsum 
     288   !!  
     289      CALL wrk_dealloc(jpi,jpj,jpk,   ztmp3d )  
     290      CALL wrk_dealloc(jpi,jpj,       zde3t  )  
     291      CALL wrk_dealloc(jpi,jpj,       zssh0  )  
    318292   END SUBROUTINE iscpl_cons 
    319293 
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90

    r5790 r5802  
    5959      CHARACTER(20) :: cfile 
    6060 
    61       CALL wrk_alloc(jpi,jpj,jpk, ztmask_b, zumask_b, zvmask_b) ! need this new variables for evoving isf cavity 
    62       CALL wrk_alloc(jpi,jpj,jpk, ze3t_b  , ze3u_b  , ze3v_b  ) ! need of this variables for interpolation) 
     61      CALL wrk_alloc(jpi,jpj,jpk, ztmask_b, zumask_b, zvmask_b) ! mask before 
     62      CALL wrk_alloc(jpi,jpj,jpk, ze3t_b  , ze3u_b  , ze3v_b  ) ! e3   before 
    6363      CALL wrk_alloc(jpi,jpj,jpk, zdepw_b ) 
    6464      CALL wrk_alloc(jpi,jpj,     zsmask_b                    ) 
     
    105105      CALL wrk_dealloc(jpi,jpj,jpk, zdepw_b                    ) 
    106106      CALL wrk_dealloc(jpi,jpj,     zsmask_b                   ) 
    107   
    108       neuler = 0              ! Euler restart (neuler=0) 
    109  
    110       tsb    (:,:,:,:) = tsn    (:,:,:,:)  ! all before fields set to now values 
    111       ub     (:,:,:)   = un     (:,:,:) 
    112       vb     (:,:,:)   = vn     (:,:,:) 
    113       rotb   (:,:,:)   = rotn   (:,:,:) 
    114       hdivb  (:,:,:)   = hdivn  (:,:,:) 
    115       fse3t_b(:,:,:)   = fse3t_n(:,:,:) 
    116       sshb   (:,:)     = sshn   (:,:) 
     107 
     108      !! next step is an euler time step 
     109      neuler = 0 
     110      !! set _b and _n variables equal 
     111      tsb (:,:,:,:) = tsn (:,:,:,:) 
     112      ub  (:,:,:  ) = un  (:,:,:  ) 
     113      vb  (:,:,:  ) = vn  (:,:,:  ) 
     114      sshb(:,:    ) = sshn(:,:) 
     115      !! set _b and _n vertical scale factor equal 
     116      fse3t_b (:,:,:) = fse3t_n (:,:,:) 
     117      fse3u_b (:,:,:) = fse3u_n (:,:,:) 
     118      fse3v_b (:,:,:) = fse3v_n (:,:,:) 
     119      IF ( lk_vvl ) THEN 
     120         fse3uw_b(:,:,:) = fse3uw_n(:,:,:) 
     121         fse3vw_b(:,:,:) = fse3vw_n(:,:,:) 
     122         fsdept_b(:,:,:) = fsdept_n(:,:,:) 
     123         fsdepw_b(:,:,:) = fsdepw_n(:,:,:) 
     124         hu_b (:,:) = hu(:,:) 
     125         hv_b (:,:) = hv(:,:) 
     126         hur_b(:,:) = hur(:,:) 
     127         hvr_b(:,:) = hvr(:,:) 
     128      END IF 
    117129      ! 
    118130   END SUBROUTINE iscpl_stp 
     
    166178      END DO 
    167179            
    168       ! compute ssh      
    169       zsmask0(:,:) = psmask_b(:,:) 
    170       zsmask1(:,:) = psmask_b(:,:) 
    171        
    172180      ! compute new ssh if we open a full water column (average of the closest neigbourgs)   
    173181      sshb (:,:)=sshn(:,:) 
    174182      zssh0(:,:)=sshn(:,:) 
     183      zsmask0(:,:) = psmask_b(:,:) 
     184      zsmask1(:,:) = psmask_b(:,:) 
    175185      DO iz = 1,10    ! need to be tuned (configuration dependent) 
    176186         zdsmask(:,:) = ssmask(:,:)-zsmask0(:,:) 
     
    198208 
    199209!============================================================================= 
    200 ! WARNING we cannot used glob_sum for before time variable (because glob_sum use tmask_i). Need to mask the halo and glob_sum_full 
    201       ztmp3d(:,:,:) = 0.0 
    202       ztmp3d(2:jpi-1,2:jpj-1,:) = fse3t_n(2:jpi-1,2:jpj-1,:) * tmask(2:jpi-1,2:jpj-1,:) & 
    203       &                         - pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:) 
    204       zsum = glob_sum_full(ztmp3d)  
    205       IF (lwp) PRINT *, 'total volume correction 00 = ',zsum 
    206210      IF ( lk_vvl ) THEN 
    207211! compute fse3t_n 
    208212         DO jk = 1,jpk 
    209             fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1 + sshn(:,:) / ( ht_0(:,:) + 1. - ssmask(:,:) ) * tmask(:,:,jk) ) 
    210          END DO 
    211          ztmp3d(:,:,:) = 0.0 
    212          ztmp3d(2:jpi-1,2:jpj-1,:) = fse3t_n(2:jpi-1,2:jpj-1,:) * tmask(2:jpi-1,2:jpj-1,:) & 
    213          &                         - pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:) 
    214          zsum = glob_sum_full(ztmp3d)  
    215          IF (lwp) PRINT *, 'total volume correction 01 = ',zsum 
     213            fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:) / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) ) 
     214         END DO 
    216215 
    217216! compute fse3u/v ... (call interpolation vvl) 
     
    269268 
    270269      END IF 
    271       ztmp3d(:,:,:) = 0.0 
    272       ztmp3d(2:jpi-1,2:jpj-1,:) = fse3t_n(2:jpi-1,2:jpj-1,:) * tmask(2:jpi-1,2:jpj-1,:) & 
    273       &                         - pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:) 
    274       zsum = glob_sum_full(ztmp3d)  
    275       IF (lwp) PRINT *, 'total volume correction 02 = ',zsum 
    276270 
    277271!============================================================================= 
     
    289283         END DO 
    290284      END DO 
     285 
    291286! compute new velocity if we close a cell (check barotropic velocity and change velocity over the water column) 
    292287! compute barotropic velocity now and after  
     
    299294      ztrp(:,:,:) = vn(:,:,:)*fse3v_n(:,:,:);  
    300295      zbvn(:,:)   = SUM(ztrp,DIM=3) 
     296 
    301297! Already know ???????? 
    302298      hu1=0.0_wp ; 
     
    391387      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:) 
    392388 
    393       ! CHECK heat 
    394       zsumn = glob_sum( tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) * zwmaskb(:,:,:) * zwmaskn(:,:,:)) 
    395       zsumb = glob_sum( tsb(:,:,:,jp_tem) *  pe3t_b(:,:,:) * zwmaskb(:,:,:) * zwmaskn(:,:,:)) 
    396       IF (lwp) PRINT *, 'CHECK tsn = ',zsumn, zsumb 
    397  
    398       ! compute new T/S (interpolation) because of vvl 
    399       IF ( lk_vvl .AND. .FALSE. ) THEN 
    400       !IF ( lk_vvl ) THEN 
     389      ! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask 
     390      IF ( lk_vvl ) THEN 
    401391         DO jk = 2,jpk-1 
    402392            DO jj = 1,jpj 
     
    432422      END IF 
    433423 
    434       ! CHECK heat 
    435       ztmp3d(:,:,:) = 0.0 
    436       ztmp3d(2:jpi-1,2:jpj-1,:) = fse3t_n(2:jpi-1,2:jpj-1,:) * tmask(2:jpi-1,2:jpj-1,:) & 
    437       &                         - pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:) 
    438       zsum = glob_sum_full(ztmp3d)  
    439       IF (lwp) PRINT *, 'total volume correction 03 = ',zsum 
    440  
     424      ! Special value for closed pool and set the mask to 0 to run 
    441425      WHERE (tmask(:,:,:) == 1.0 .AND. tsn(:,:,:,2) == 0._wp)  
    442426         tsn(:,:,:,2)=  -99._wp 
     
    453437      ! compute new tn and sn if we close cell  
    454438      ! nothing to do 
    455  
    456       !! next step is an euler time step 
    457       neuler = 0  
    458       !! set _b and _n variables equal 
    459       tsb (:,:,:,:)=tsn (:,:,:,:) 
    460       ub  (:,:,:  )=un  (:,:,:  ) 
    461       vb  (:,:,:  )=vn  (:,:,:  ) 
    462       sshb(:,:    )=sshn(:,:) 
    463       !! set _b and _n vertical scale factor equal 
    464       fse3t_b (:,:,:)=fse3t_n (:,:,:) 
    465       fse3u_b (:,:,:)=fse3u_n (:,:,:) 
    466       fse3v_b (:,:,:)=fse3v_n (:,:,:) 
    467       IF ( lk_vvl ) THEN 
    468          fse3uw_b(:,:,:) = fse3uw_n(:,:,:) 
    469          fse3vw_b(:,:,:) = fse3vw_n(:,:,:) 
    470          fsdept_b(:,:,:) = fsdept_n(:,:,:) 
    471          fsdepw_b(:,:,:) = fsdepw_n(:,:,:) 
    472          hu_b(:,:) = hu(:,:) 
    473          hv_b(:,:) = hv(:,:) 
    474          hur_b(:,:) = hur(:,:) 
    475          hvr_b(:,:) = hvr(:,:) 
    476       END IF 
    477       !!  
     439      !  
     440      ! deallocation tmp arrays 
    478441      CALL wrk_dealloc(jpi,jpj,jpk,2, zts0                                   ) 
    479442      CALL wrk_dealloc(jpi,jpj,jpk,   ztmask0, ztmask1 , ztrp                )  
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/lib_fortran.F90

    r5779 r5802  
    172172      !!----------------------------------------------------------------------- 
    173173      ! 
    174       glob_sum_full_2d = SUM( ptab(:,:) ) 
     174      glob_sum_full_2d = SUM( ptab(:,:) * tmask_h(:,:) ) 
    175175      IF( lk_mpp )   CALL mpp_sum( glob_sum_full_2d ) 
    176176      ! 
     
    194194      glob_sum_full_3d = 0.e0 
    195195      DO jk = 1, ijpk 
    196          glob_sum_full_3d = glob_sum_full_3d + SUM( ptab(:,:,jk) ) 
     196         glob_sum_full_3d = glob_sum_full_3d + SUM( ptab(:,:,jk) * tmask_h(:,:) ) 
    197197      END DO 
    198198      IF( lk_mpp )   CALL mpp_sum( glob_sum_full_3d ) 
     
    376376      DO jj = 1, jpj 
    377377         DO ji =1, jpi 
    378          ztmp =  ptab(ji,jj) 
     378         ztmp =  ptab(ji,jj) * tmask_h(ji,jj)  
    379379         CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 
    380380         END DO 
     
    407407         DO jj = 1, jpj 
    408408            DO ji =1, jpi 
    409             ztmp =  ptab(ji,jj,jk) 
     409            ztmp =  ptab(ji,jj,jk) * tmask_h(ji,jj) 
    410410            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 
    411411            END DO 
Note: See TracChangeset for help on using the changeset viewer.