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 5820 for branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM – NEMO

Ignore:
Timestamp:
2015-10-21T18:01:58+02:00 (8 years ago)
Author:
mathiot
Message:

ice sheet coupling: remove some print, fix pb in diahsb if ssmask is modified, rm corner extrapolation + some bug fix in conservation option

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

Legend:

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

    r5802 r5820  
    4646   REAL(wp) ::   frc_wn_t, frc_wn_s    ! global forcing trends 
    4747   ! 
    48    REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   surf          , ssh_ini          ! 
     48   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   surf  
     49   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   surf_ini      , ssh_ini          ! 
    4950   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   ssh_hc_loc_ini, ssh_sc_loc_ini   ! 
    5051   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   hc_loc_ini, sc_loc_ini, e3t_ini  ! 
     
    144145 
    145146      ! volume variation (calculated with ssh) 
    146       zdiff_v1 = glob_sum( surf(:,:) * ( sshn(:,:) - ssh_ini(:,:) ) ) 
     147      zdiff_v1 = glob_sum_full( surf(:,:) * sshn(:,:) - surf_ini(:,:) * ssh_ini(:,:) ) 
    147148 
    148149      ! heat & salt content variation (associated with ssh) 
     
    165166      DO jk = 1, jpkm1 
    166167         ! volume variation (calculated with scale factors) 
    167          zdiff_v2 = zdiff_v2 + glob_sum( surf(:,:) * ( tmask(:,:,jk) & 
    168             &                           * fse3t_n(:,:,jk) - e3t_ini(:,:,jk) ) ) 
    169          ! heat content variation 
    170          zdiff_hc = zdiff_hc + glob_sum(  surf(:,:) * ( tmask(:,:,jk) &  
    171             &                           * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) - hc_loc_ini(:,:,jk) ) ) 
     168         zdiff_v2 = zdiff_v2 + glob_sum_full( surf    (:,:) * fse3t_n  (:,:,jk) * tmask(:,:,jk) & 
     169              &                             - surf_ini(:,:) *   e3t_ini(:,:,jk) ) 
     170         ! heat content variation  
     171         zdiff_hc = zdiff_hc + glob_sum_full( surf    (:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) & 
     172              &                             - surf_ini(:,:) * hc_loc_ini(:,:,jk) ) 
    172173         ! salt content variation 
    173          zdiff_sc = zdiff_sc + glob_sum(  surf(:,:) * ( tmask(:,:,jk)   & 
    174             &                           * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) - sc_loc_ini(:,:,jk) ) ) 
     174         zdiff_sc = zdiff_sc + glob_sum_full( surf    (:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) &  
     175              &                             - surf_ini(:,:) * sc_loc_ini(:,:,jk) ) 
    175176      ENDDO 
    176177 
     
    192193      zvol_tot = 0._wp                    ! total ocean volume (calculated with scale factors) 
    193194      DO jk = 1, jpkm1 
    194          zvol_tot  = zvol_tot + glob_sum( surf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) ) 
     195         zvol_tot  = zvol_tot + glob_sum_full( surf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) ) 
    195196      END DO 
    196197 
     
    268269              CALL iom_get( numror, 'frc_wn_s', frc_wn_s ) 
    269270           ENDIF 
     271           CALL iom_get( numror, jpdom_autoglo, 'surf_ini', surf_ini ) 
    270272           CALL iom_get( numror, jpdom_autoglo, 'ssh_ini', ssh_ini ) 
    271273           CALL iom_get( numror, jpdom_autoglo, 'e3t_ini', e3t_ini ) 
     
    280282          IF(lwp) WRITE(numout,*) ' dia_hsb at initial state ' 
    281283          IF(lwp) WRITE(numout,*) '~~~~~~~' 
     284          surf_ini(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:)             ! initial ocean surface 
    282285          ssh_ini(:,:) = sshn(:,:)                                       ! initial ssh 
    283286          DO jk = 1, jpk 
     
    320323           CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_s', frc_wn_s ) 
    321324        ENDIF 
     325        CALL iom_rstput( kt, nitrst, numrow, 'surf_ini', surf_ini ) 
    322326        CALL iom_rstput( kt, nitrst, numrow, 'ssh_ini', ssh_ini ) 
    323327        CALL iom_rstput( kt, nitrst, numrow, 'e3t_ini', e3t_ini ) 
     
    387391      ! 1 - Allocate memory ! 
    388392      ! ------------------- ! 
    389       ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), & 
    390          &      e3t_ini(jpi,jpj,jpk), surf(jpi,jpj),  ssh_ini(jpi,jpj), STAT=ierror ) 
     393      ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), surf_ini(jpi,jpj), & 
     394         &      e3t_ini(jpi,jpj,jpk), surf(jpi,jpj),  ssh_ini(jpi,jpj), STAT=ierror  ) 
    391395      IF( ierror > 0 ) THEN 
    392396         CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN 
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90

    r5802 r5820  
    6262      REAL(wp):: summsk, zsum, zsum1, zarea, zsumn, zsumb 
    6363      REAL(wp):: r1_tiscpl 
    64       REAL(wp):: zjip1_ratio, zjim1_ratio, zjjp1_ratio, zjjm1_ratio 
     64      REAL(wp):: zjip1_ratio  , zjim1_ratio  , zjjp1_ratio  , zjjm1_ratio 
    6565      !! 
    6666      REAL(wp), DIMENSION(:,:    ), POINTER :: zde3t, zdtem, zdsal 
     
    8383      pts_flx (:,:,:,:) = 0.0_wp 
    8484 
    85       ! mask tsn and tsb (should be useless) 
     85      ! mask tsn and tsb  
    8686      tsb(:,:,:,jp_tem)=tsb(:,:,:,jp_tem)*ptmask_b(:,:,:); tsn(:,:,:,jp_tem)=tsn(:,:,:,jp_tem)*tmask(:,:,:); 
    8787      tsb(:,:,:,jp_sal)=tsb(:,:,:,jp_sal)*ptmask_b(:,:,:); tsn(:,:,:,jp_sal)=tsn(:,:,:,jp_sal)*tmask(:,:,:); 
     
    9898               IF (tmask_h(ji,jj) == 1._wp) THEN 
    9999 
    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)   & 
     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)   & 
    105105                            - 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) 
     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) 
    109109                
    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) 
     110                  ! shh changes 
     111                  IF ( ptmask_b(ji,jj,jk) == 1._wp .OR. tmask(ji,jj,jk) == 1._wp ) THEN  
     112                     zde3t(ji,jj) = zde3t(ji,jj) + zssh0(ji,jj) ! zssh0 = 0 if vvl 
     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 
    160147 
    161148                        ! set to 0 the cell we distributed over neigbourg cells 
     
    163150                        pts_flx (ji,jj,jk,jp_sal) = 0._wp 
    164151                        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. 
     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                           WRITE(numout,*) 'W A R N I N G iscpl: no adjacent cell on the vertical and horizontal' 
     168                           WRITE(numout,*) '                     ',mig(ji),' ',mjg(jj),' ',jk 
     169                           WRITE(numout,*) '                     ',ji,' ',jj,' ',jk,' ',narea 
     170                           WRITE(numout,*) ' we are now looking for the closest wet cell on the horizontal ' 
     171                        ! We deal with these points later. 
     172                        END IF 
    172173                     END IF 
    173174                  END IF 
    174175               END IF 
    175                END IF 
    176176            END DO 
    177177         END DO 
     
    181181      CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1.) 
    182182      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 
    190183 
    191184      ! if no neighbour wet cell in case of 2close a cell", need to find the nearest wet point  
     
    198191         DO ji = 2,jpi-1 
    199192            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 
     193               IF (     ptmask_b(ji,jj,jk) == 1._wp .AND. tmask(ji,jj,jk+1)  == 0._wp .AND. tmask_h(ji,jj) == 1._wp  & 
     194                  .AND. SUM(tmask(ji-1:ji+1,jj,jk)) + SUM(tmask(ji,jj-1:jj+1,jk)) == 0._wp) THEN 
    203195                  vnpts(narea) = vnpts(narea) + 1  
    204196               END IF 
     
    216208      ALLOCATE(ixpts(npts), iypts(npts), izpts(npts), zcorr_vol(npts), zcorr_sal(npts), zcorr_tem(npts), zlon(npts), zlat(npts)) 
    217209      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 
     210      zcorr_vol(:) = -1.0e20 
     211      zcorr_sal(:) = -1.0e20 
     212      zcorr_tem(:) = -1.0e20 
    221213 
    222214      ! fill new variable 
     
    225217         DO ji = 2,jpi-1 
    226218            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 
     219               IF (     ptmask_b(ji,jj,jk) == 1._wp .AND. tmask(ji,jj,jk+1)  == 0._wp .AND. tmask_h(ji,jj) == 1._wp  & 
     220                  .AND. SUM(tmask(ji-1:ji+1,jj,jk)) + SUM(tmask(ji,jj-1:jj+1,jk)) == 0._wp) THEN 
    230221                  jpts = jpts + 1  ! positioning in the vnpts vector for the area narea 
    231                   PRINT *, 'corrected point ', narea, ji, jj, jk, jpts 
    232222                  ixpts(jpts) = ji           ; iypts(jpts) = jj ; izpts(jpts) = jk 
    233223                  zlon (jpts) = glamt(ji,jj) ; zlat (jpts) = gphit(ji,jj) 
     
    235225                  zcorr_sal(jpts) = pts_flx (ji,jj,jk,jp_sal) 
    236226                  zcorr_tem(jpts) = pts_flx (ji,jj,jk,jp_tem) 
     227 
    237228                  ! set flx to 0 (safer) 
    238229                  pvol_flx(ji,jj,jk       ) = 0.0_wp 
    239230                  pts_flx (ji,jj,jk,jp_sal) = 0.0_wp 
    240231                  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 
    242232               END IF 
    243233            END DO 
     
    246236 
    247237      ! build array of total problematic point on each cpu (share to each cpu) 
     238      ! point coordinates 
    248239      CALL mpp_max(zlat ,npts) 
    249240      CALL mpp_max(zlon ,npts) 
    250241      CALL mpp_max(izpts,npts) 
    251242 
     243      ! correction values  
     244      CALL mpp_max(zcorr_vol,npts) 
     245      CALL mpp_max(zcorr_sal,npts) 
     246      CALL mpp_max(zcorr_tem,npts) 
     247 
    252248      ! put correction term in the closest cell           
    253       PRINT *, 'corrected point1 ', narea, zlon, zlat, izpts 
    254249      DO jpts = 1,npts 
    255250         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) 
    257251         DO jj = mj0(iypts(jpts)),mj1(iypts(jpts)) 
    258252            DO ji = mi0(ixpts(jpts)),mi1(ixpts(jpts)) 
    259253               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 
     254 
     255               IF (tmask_h(ji,jj) == 1._wp) THEN 
     256                  ! correct the vol_flx in the closest cell 
     257                  pvol_flx(ji,jj,jk)        =  pvol_flx(ji,jj,jk       ) + zcorr_vol(jpts) 
     258                  pts_flx (ji,jj,jk,jp_sal) =  pts_flx (ji,jj,jk,jp_sal) + zcorr_sal(jpts) 
     259                  pts_flx (ji,jj,jk,jp_tem) =  pts_flx (ji,jj,jk,jp_tem) + zcorr_tem(jpts) 
     260 
     261                  ! set correction to 0 
     262                  zcorr_vol(jpts) = 0.0_wp 
     263                  zcorr_sal(jpts) = 0.0_wp 
     264                  zcorr_tem(jpts) = 0.0_wp 
     265               END IF 
     266            END DO 
     267         END DO 
     268      END DO 
     269 
    266270      ! deallocate variables  
    267271      DEALLOCATE(vnpts) 
     
    273277      pts_flx (:,:,:,jp_tem) = pts_flx (:,:,:,jp_tem) * tmask(:,:,:) 
    274278 
    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    !!  
     279      ! compute sum over the halo and set it to 0. 
     280      CALL lbc_sum(pvol_flx(:,:,:       ),'T',1._wp) 
     281      CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1._wp) 
     282      CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1._wp) 
     283 
     284      ! deallocate variables 
    289285      CALL wrk_dealloc(jpi,jpj,jpk,   ztmp3d )  
    290286      CALL wrk_dealloc(jpi,jpj,       zde3t  )  
  • branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90

    r5802 r5820  
    9797         CALL iom_rstput( 0, 0, inum0, 'tem_cor', htsc_iscpl(:,:,:,jp_tem) ) 
    9898         CALL iom_rstput( 0, 0, inum0, 'sal_cor', htsc_iscpl(:,:,:,jp_sal) ) 
    99          CALL iom_rstput( 0, 0, inum0, 'e3t_diff', (e3t_0(:,:,:)*tmask(:,:,:) - ze3t_b(:,:,:)*ztmask_b(:,:,:))*tmask(:,:,:) ) 
    10099         CALL iom_close ( inum0 ) 
    101100      END IF 
     
    191190               summsk=(zsmask0(jip1,jj)+zsmask0(jim1,jj)+zsmask0(ji,jjp1)+zsmask0(ji,jjm1)) 
    192191               IF (zdsmask(ji,jj)==1 .AND. summsk .NE. 0) THEN 
    193                   PRINT *, 'add ssh to 1 cell',ji,jj,narea 
    194                   sshn(ji,jj)=(zssh0(jip1,jj)*zsmask0(jip1,jj)     & 
    195                   &           +zssh0(jim1,jj)*zsmask0(jim1,jj)     & 
    196                   &           +zssh0(ji,jjp1)*zsmask0(ji,jjp1)     & 
    197                   &           +zssh0(ji,jjm1)*zsmask0(ji,jjm1))/summsk 
     192                  sshn(ji,jj)=( zssh0(jip1,jj)*zsmask0(jip1,jj)     & 
     193                  &           + zssh0(jim1,jj)*zsmask0(jim1,jj)     & 
     194                  &           + zssh0(ji,jjp1)*zsmask0(ji,jjp1)     & 
     195                  &           + zssh0(ji,jjm1)*zsmask0(ji,jjm1))/summsk 
    198196                  zsmask1(ji,jj)=1 
    199197               END IF 
     
    209207!============================================================================= 
    210208      IF ( lk_vvl ) THEN 
    211 ! compute fse3t_n 
    212          DO jk = 1,jpk 
    213             fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:) / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) ) 
    214          END DO 
    215  
    216 ! compute fse3u/v ... (call interpolation vvl) 
    217       ! Reconstruction of all vertical scale factors at now and before time steps 
     209      ! Reconstruction of all vertical scale factors at now time steps 
    218210      ! ============================================================================= 
    219211      ! Horizontal scale factor interpolations 
    220212      ! -------------------------------------- 
     213         DO jk = 1,jpk 
     214            DO jj=1,jpj 
     215               DO ji=1,jpi 
     216                  IF (tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp) THEN 
     217                     fse3t_n(ji,jj,jk) = e3t_0(ji,jj,jk) * ( 1._wp + sshn(ji,jj) / ( ht_0(ji,jj) + 1._wp - ssmask(ji,jj) ) * tmask(ji,jj,jk) ) 
     218                  ENDIF 
     219               END DO 
     220            END DO 
     221         END DO 
     222 
    221223         CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 
    222224         CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 
     
    295297      zbvn(:,:)   = SUM(ztrp,DIM=3) 
    296298 
    297 ! Already know ???????? 
     299      ! new water column 
    298300      hu1=0.0_wp ; 
    299301      hv1=0.0_wp ; 
     
    302304        hv1(:,:) = hv1(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk) 
    303305      END DO 
     306       
     307      ! compute correction       
    304308      zucorr = 0._wp 
    305309      zvcorr = 0._wp 
     
    314318         END DO 
    315319      END DO  
     320       
     321      ! update velocity 
    316322      DO jk  = 1,jpk 
    317323         un(:,:,jk)=(un(:,:,jk) - zucorr(:,:))*umask(:,:,jk) 
     
    333339                      jip1=ji+1; jim1=ji-1; 
    334340                      jjp1=jj+1; jjm1=jj-1; 
    335                       summsk=(ztmask0(jip1,jj,jk)+ztmask0(jim1,jj,jk)+ztmask0(ji,jjp1,jk)+ztmask0(ji,jjm1,jk)) 
    336                       !! horizontal extrapolation 
     341                      summsk= (ztmask0(jip1,jj  ,jk)+ztmask0(jim1,jj  ,jk)+ztmask0(ji  ,jjp1,jk)+ztmask0(ji  ,jjm1,jk)) 
     342                      !! horizontal basic extrapolation 
    337343                      IF (zdmask(ji,jj)==1 .AND. summsk .NE. 0) THEN 
    338                          tsn(ji,jj,jk,1)=( zts0(jip1,jj,jk,1)*ztmask0(jip1,jj,jk)    & 
    339                          &                +zts0(jim1,jj,jk,1)*ztmask0(jim1,jj,jk)    & 
    340                          &                +zts0(ji,jjp1,jk,1)*ztmask0(ji,jjp1,jk)    & 
    341                          &                +zts0(ji,jjm1,jk,1)*ztmask0(ji,jjm1,jk))/summsk 
    342                          tsn(ji,jj,jk,2)=( zts0(jip1,jj,jk,2)*ztmask0(jip1,jj,jk)    & 
    343                          &                +zts0(jim1,jj,jk,2)*ztmask0(jim1,jj,jk)    & 
    344                          &                +zts0(ji,jjp1,jk,2)*ztmask0(ji,jjp1,jk)    & 
    345                          &                +zts0(ji,jjm1,jk,2)*ztmask0(ji,jjm1,jk))/summsk 
     344                         tsn(ji,jj,jk,1)=( zts0(jip1,jj  ,jk,1)*ztmask0(jip1,jj  ,jk) & 
     345                         &                +zts0(jim1,jj  ,jk,1)*ztmask0(jim1,jj  ,jk) & 
     346                         &                +zts0(ji  ,jjp1,jk,1)*ztmask0(ji  ,jjp1,jk) & 
     347                         &                +zts0(ji  ,jjm1,jk,1)*ztmask0(ji  ,jjm1,jk) ) / summsk 
     348                         tsn(ji,jj,jk,2)=( zts0(jip1,jj  ,jk,2)*ztmask0(jip1,jj  ,jk) & 
     349                         &                +zts0(jim1,jj  ,jk,2)*ztmask0(jim1,jj  ,jk) & 
     350                         &                +zts0(ji  ,jjp1,jk,2)*ztmask0(ji  ,jjp1,jk) & 
     351                         &                +zts0(ji  ,jjm1,jk,2)*ztmask0(ji  ,jjm1,jk) ) / summsk 
    346352                         ztmask1(ji,jj,jk)=1 
    347353                      END IF 
     
    358364                         END IF 
    359365                      END IF 
    360                       !! horizontal corner extrapolation if horizontal and vertical failed 
    361                       IF (zdmask(ji,jj)==1 .AND. summsk==0) THEN 
    362                          jip1=ji+1; jim1=ji-1; 
    363                          jjp1=jj+1; jjm1=jj-1; 
    364                          summsk=(ztmask0(jip1,jjp1,jk)+ztmask0(jim1,jjm1,jk)+ztmask0(jim1,jjp1,jk)+ztmask0(jip1,jjm1,jk)) 
    365                          IF (zdmask(ji,jj)==1 .AND. summsk .NE. 0) THEN 
    366                             tsn(ji,jj,jk,1)=( zts0(jip1,jjp1,jk,1)*ztmask0(jip1,jjp1,jk)     & 
    367                             &                +zts0(jim1,jjm1,jk,1)*ztmask0(jim1,jjm1,jk)     & 
    368                             &                +zts0(jim1,jjp1,jk,1)*ztmask0(jim1,jjp1,jk)     & 
    369                             &                +zts0(jip1,jjm1,jk,1)*ztmask0(jip1,jjm1,jk))/summsk 
    370                             tsn(ji,jj,jk,2)=( zts0(jip1,jjp1,jk,2)*ztmask0(jip1,jjp1,jk)     & 
    371                             &                +zts0(jim1,jjm1,jk,2)*ztmask0(jim1,jjm1,jk)     & 
    372                             &                +zts0(jim1,jjp1,jk,2)*ztmask0(jim1,jjp1,jk)     & 
    373                             &                +zts0(jip1,jjm1,jk,2)*ztmask0(jip1,jjm1,jk))/summsk 
    374                             ztmask1(ji,jj,jk)=1 
    375                          END IF 
    376                       END IF 
    377366                  END DO 
    378367              END DO 
     
    392381            DO jj = 1,jpj 
    393382               DO ji = 1,jpi 
    394                   IF (zwmaskn(ji,jj,jk) * zwmaskb(ji,jj,jk) == 1.0_wp ) THEN 
     383                  IF (zwmaskn(ji,jj,jk) * zwmaskb(ji,jj,jk) == 1.0_wp .AND. (tmask(ji,jj,1)==0 .OR. ptmask_b(ji,jj,1)==0) ) THEN 
    395384                     !compute weight 
    396385                     zdzp1 = MAX(0._wp,fsdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk+1)) 
     
    398387                     zdzm1 = MAX(0._wp,pdepw_b(ji,jj,jk  ) - fsdepw_n(ji,jj,jk  )) 
    399388                     IF (zdz .LT. 0.0_wp) THEN  
    400                         PRINT *, 'ERROR dz n ', ji,jj,jk,zdz,fsdepw_n(ji,jj,jk+1),fsdepw_n(ji,jj,jk),fsdepw_n(ji,jj,jk-1) 
    401                         PRINT *, 'ERROR dz n             = ',fse3t_n (ji,jj,jk+1),fse3t_n (ji,jj,jk),fse3t_n (ji,jj,jk-1), sshn(ji,jj) 
    402                         PRINT *, 'ERROR dz b             = ',pdepw_b(ji,jj,jk+1),pdepw_b(ji,jj,jk),pdepw_b(ji,jj,jk-1) 
    403                         PRINT *, 'ERROR dz b             = ',fse3t_b (ji,jj,jk+1),fse3t_b (ji,jj,jk),fse3t_b (ji,jj,jk-1), sshb(ji,jj) 
    404                         PRINT *, 'ERROR dz 0             = ',  e3t_0 (ji,jj,jk+1),  e3t_0 (ji,jj,jk),  e3t_0 (ji,jj,jk-1) 
    405                         PRINT *, 'ERROR dz n             = ',  tmask (ji,jj,jk+1),  tmask (ji,jj,jk),  tmask (ji,jj,jk-1) 
    406                         PRINT *, 'ERROR dz n             = ', zwmaskn(ji,jj,jk+1), zwmaskn(ji,jj,jk), zwmaskn(ji,jj,jk-1) 
    407                         PRINT *, 'ERROR dz b             = ', ptmask_b(ji,jj,jk+1), ptmask_b(ji,jj,jk), ptmask_b(ji,jj,jk-1) 
    408                         PRINT *, 'ERROR dz b             = ', zwmaskb(ji,jj,jk+1), zwmaskb(ji,jj,jk), zwmaskb(ji,jj,jk-1) 
    409                         PRINT *, 'ERROR dz b             = ', gdepw_0(ji,jj,jk+1), gdepw_0(ji,jj,jk), gdepw_0(ji,jj,jk-1) 
    410                         STOP 
     389                        WRITE(numout,*) 'ERROR dz n ', ji,jj,jk,zdz,fsdepw_n(ji,jj,jk+1),fsdepw_n(ji,jj,jk),fsdepw_n(ji,jj,jk-1) 
     390                        WRITE(numout,*) 'ERROR dz n             = ',fse3t_n (ji,jj,jk+1),fse3t_n (ji,jj,jk),fse3t_n (ji,jj,jk-1), sshn(ji,jj) 
     391                        WRITE(numout,*) 'ERROR dz b             = ',pdepw_b(ji,jj,jk+1),pdepw_b(ji,jj,jk),pdepw_b(ji,jj,jk-1) 
     392                        WRITE(numout,*) 'ERROR dz b             = ',fse3t_b (ji,jj,jk+1),fse3t_b (ji,jj,jk),fse3t_b (ji,jj,jk-1), sshb(ji,jj) 
     393                        WRITE(numout,*) 'ERROR dz 0             = ',  e3t_0 (ji,jj,jk+1),  e3t_0 (ji,jj,jk),  e3t_0 (ji,jj,jk-1) 
     394                        WRITE(numout,*) 'ERROR dz n             = ',  tmask (ji,jj,jk+1),  tmask (ji,jj,jk),  tmask (ji,jj,jk-1) 
     395                        WRITE(numout,*) 'ERROR dz n             = ', zwmaskn(ji,jj,jk+1), zwmaskn(ji,jj,jk), zwmaskn(ji,jj,jk-1) 
     396                        WRITE(numout,*) 'ERROR dz b             = ', ptmask_b(ji,jj,jk+1), ptmask_b(ji,jj,jk), ptmask_b(ji,jj,jk-1) 
     397                        WRITE(numout,*) 'ERROR dz b             = ', zwmaskb(ji,jj,jk+1), zwmaskb(ji,jj,jk), zwmaskb(ji,jj,jk-1) 
     398                        WRITE(numout,*) 'ERROR dz b             = ', gdepw_0(ji,jj,jk+1), gdepw_0(ji,jj,jk), gdepw_0(ji,jj,jk-1) 
     399                        CALL ctl_stop( 'STOP', 'rst_iscpl : unable to compute the interpolation' ) 
    411400                     END IF 
    412401                     tsn(ji,jj,jk,jp_tem) = ( zdzp1*tsb(ji,jj,jk+1,jp_tem)      & 
Note: See TracChangeset for help on using the changeset viewer.