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 11730 for NEMO/releases/release-4.0.1/src/ICE/icedyn_adv_pra.F90 – NEMO

Ignore:
Timestamp:
2019-10-18T17:47:20+02:00 (5 years ago)
Author:
clem
Message:

last commit: avoid division by 0 occuring on some rare occasions (CMMC and ECMWF issues). It might be machine dependant.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • NEMO/releases/release-4.0.1/src/ICE/icedyn_adv_pra.F90

    r11632 r11730  
    3939   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxsal, sysal, sxxsal, syysal, sxysal   ! ice salinity 
    4040   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxage, syage, sxxage, syyage, sxyage   ! ice age 
    41    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxopw, syopw, sxxopw, syyopw, sxyopw   ! open water in sea ice 
    4241   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sxc0 , syc0 , sxxc0 , syyc0 , sxyc0    ! snow layers heat content 
    4342   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sxe  , sye  , sxxe  , syye  , sxye     ! ice layers heat content 
    … …  
    8180      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
    8281      ! 
    83       INTEGER  ::   jk, jl, jt              ! dummy loop indices 
     82      INTEGER  ::   ji,jj, jk, jl, jt       ! dummy loop indices 
    8483      INTEGER  ::   icycle                  ! number of sub-timestep for the advection 
    8584      REAL(wp) ::   zdt                     !   -      - 
    8685      REAL(wp), DIMENSION(1)                  ::   zcflprv, zcflnow   ! for global communication 
     86      REAL(wp), DIMENSION(jpi,jpj)            ::   zati1, zati2 
     87      REAL(wp), DIMENSION(jpi,jpj)            ::   zudy, zvdx 
    8788      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zarea 
    88       REAL(wp), DIMENSION(jpi,jpj,1)          ::   z0opw 
    8989      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ice, z0snw, z0ai, z0smi, z0oi 
    9090      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp 
    … …  
    109109      zdt = rdt_ice / REAL(icycle) 
    110110       
    111       !------------------------- 
    112       ! transported fields                                         
    113       !------------------------- 
    114       z0opw(:,:,1) = pato_i(:,:) * e1e2t(:,:)              ! Open water area  
    115       DO jl = 1, jpl 
    116          zarea(:,:,jl) = e1e2t(:,:) 
    117          z0snw(:,:,jl) = pv_s (:,:,jl) * e1e2t(:,:)        ! Snow volume 
    118          z0ice(:,:,jl) = pv_i (:,:,jl) * e1e2t(:,:)        ! Ice  volume 
    119          z0ai (:,:,jl) = pa_i (:,:,jl) * e1e2t(:,:)        ! Ice area 
    120          z0smi(:,:,jl) = psv_i(:,:,jl) * e1e2t(:,:)        ! Salt content 
    121          z0oi (:,:,jl) = poa_i(:,:,jl) * e1e2t(:,:)        ! Age content 
    122          DO jk = 1, nlay_s 
    123             z0es(:,:,jk,jl) = pe_s(:,:,jk,jl) * e1e2t(:,:) ! Snow heat content 
    124          END DO 
    125          DO jk = 1, nlay_i 
    126             z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content 
    127          END DO 
    128          IF ( ln_pnd_H12 ) THEN 
    129             z0ap(:,:,jl)  = pa_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond fraction 
    130             z0vp(:,:,jl)  = pv_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond volume 
    131          ENDIF 
    132       END DO 
    133  
    134       !                                                    !--------------------------------------------! 
    135       IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==! 
    136          !                                                 !--------------------------------------------! 
    137          DO jt = 1, icycle 
    138             CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0opw , sxopw , sxxopw , syopw , syyopw , sxyopw ) !--- open water 
    139             CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0opw , sxopw , sxxopw , syopw , syyopw , sxyopw ) 
    140             CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume 
    141             CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) 
    142             CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) !--- snow volume 
    143             CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) 
    144             CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity 
    145             CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) 
    146             CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) !--- ice concentration 
    147             CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) 
    148             CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) !--- ice age 
    149             CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) 
     111      ! --- transport --- ! 
     112      zudy(:,:) = pu_ice(:,:) * e2u(:,:) 
     113      zvdx(:,:) = pv_ice(:,:) * e1v(:,:) 
     114 
     115      DO jt = 1, icycle 
     116 
     117         ! record at_i before advection (for open water) 
     118         zati1(:,:) = SUM( pa_i(:,:,:), dim=3 ) 
     119          
     120         ! --- transported fields --- !                                         
     121         DO jl = 1, jpl 
     122            zarea(:,:,jl) = e1e2t(:,:) 
     123            z0snw(:,:,jl) = pv_s (:,:,jl) * e1e2t(:,:)        ! Snow volume 
     124            z0ice(:,:,jl) = pv_i (:,:,jl) * e1e2t(:,:)        ! Ice  volume 
     125            z0ai (:,:,jl) = pa_i (:,:,jl) * e1e2t(:,:)        ! Ice area 
     126            z0smi(:,:,jl) = psv_i(:,:,jl) * e1e2t(:,:)        ! Salt content 
     127            z0oi (:,:,jl) = poa_i(:,:,jl) * e1e2t(:,:)        ! Age content 
     128            DO jk = 1, nlay_s 
     129               z0es(:,:,jk,jl) = pe_s(:,:,jk,jl) * e1e2t(:,:) ! Snow heat content 
     130            END DO 
     131            DO jk = 1, nlay_i 
     132               z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content 
     133            END DO 
     134            IF ( ln_pnd_H12 ) THEN 
     135               z0ap(:,:,jl)  = pa_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond fraction 
     136               z0vp(:,:,jl)  = pv_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond volume 
     137            ENDIF 
     138         END DO 
     139         ! 
     140         !                                                                  !--------------------------------------------! 
     141         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==! 
     142            !                                                               !--------------------------------------------! 
     143            CALL adv_x( zdt , zudy , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume 
     144            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) 
     145            CALL adv_x( zdt , zudy , 1._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) !--- snow volume 
     146            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) 
     147            CALL adv_x( zdt , zudy , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity 
     148            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) 
     149            CALL adv_x( zdt , zudy , 1._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) !--- ice concentration 
     150            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) 
     151            CALL adv_x( zdt , zudy , 1._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) !--- ice age 
     152            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) 
    150153            ! 
    151             DO jk = 1, nlay_s                                                                             !--- snow heat content 
    152                CALL adv_x( zdt, pu_ice, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   & 
    153                   &                                   sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 
    154                CALL adv_y( zdt, pv_ice, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   & 
    155                   &                                   sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 
    156             END DO 
    157             DO jk = 1, nlay_i                                                                             !--- ice heat content 
    158                CALL adv_x( zdt, pu_ice, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   &  
    159                   &                                   sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 
    160                CALL adv_y( zdt, pv_ice, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   &  
    161                   &                                   sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 
     154            DO jk = 1, nlay_s                                                                           !--- snow heat content 
     155               CALL adv_x( zdt, zudy, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   & 
     156                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 
     157               CALL adv_y( zdt, zvdx, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   & 
     158                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 
     159            END DO 
     160            DO jk = 1, nlay_i                                                                           !--- ice heat content 
     161               CALL adv_x( zdt, zudy, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   &  
     162                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 
     163               CALL adv_y( zdt, zvdx, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   &  
     164                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 
    162165            END DO 
    163166            ! 
    164167            IF ( ln_pnd_H12 ) THEN 
    165                CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction 
    166                CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )  
    167                CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume 
    168                CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )  
     168               CALL adv_x( zdt , zudy , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction 
     169               CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )  
     170               CALL adv_x( zdt , zudy , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume 
     171               CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )  
    169172            ENDIF 
    170          END DO 
    171       !                                                    !--------------------------------------------! 
    172       ELSE                                                 !== even ice time step:  adv_y then adv_x  ==! 
    173          !                                                 !--------------------------------------------! 
    174          DO jt = 1, icycle 
    175             CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0opw , sxopw , sxxopw , syopw , syyopw , sxyopw ) !--- open water 
    176             CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0opw , sxopw , sxxopw , syopw , syyopw , sxyopw ) 
    177             CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume 
    178             CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) 
    179             CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) !--- snow volume 
    180             CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) 
    181             CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity 
    182             CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) 
    183             CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) !--- ice concentration 
    184             CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) 
    185             CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) !--- ice age 
    186             CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) 
    187             DO jk = 1, nlay_s                                                                             !--- snow heat content 
    188                CALL adv_y( zdt, pv_ice, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   & 
    189                   &                                   sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 
    190                CALL adv_x( zdt, pu_ice, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   & 
    191                   &                                   sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 
    192             END DO 
    193             DO jk = 1, nlay_i                                                                             !--- ice heat content 
    194                CALL adv_y( zdt, pv_ice, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   &  
    195                   &                                   sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 
    196                CALL adv_x( zdt, pu_ice, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   &  
    197                   &                                   sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 
     173            !                                                               !--------------------------------------------! 
     174         ELSE                                                               !== even ice time step:  adv_y then adv_x  ==! 
     175            !                                                               !--------------------------------------------! 
     176            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume 
     177            CALL adv_x( zdt , zudy , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) 
     178            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) !--- snow volume 
     179            CALL adv_x( zdt , zudy , 0._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) 
     180            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity 
     181            CALL adv_x( zdt , zudy , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) 
     182            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) !--- ice concentration 
     183            CALL adv_x( zdt , zudy , 0._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) 
     184            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) !--- ice age 
     185            CALL adv_x( zdt , zudy , 0._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) 
     186            DO jk = 1, nlay_s                                                                           !--- snow heat content 
     187               CALL adv_y( zdt, zvdx, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   & 
     188                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 
     189               CALL adv_x( zdt, zudy, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   & 
     190                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 
     191            END DO 
     192            DO jk = 1, nlay_i                                                                           !--- ice heat content 
     193               CALL adv_y( zdt, zvdx, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   &  
     194                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 
     195               CALL adv_x( zdt, zudy, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   &  
     196                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 
    198197            END DO 
    199198            IF ( ln_pnd_H12 ) THEN 
    200                CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction 
    201                CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 
    202                CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume 
    203                CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 
     199               CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction 
     200               CALL adv_x( zdt , zudy , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 
     201               CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume 
     202               CALL adv_x( zdt , zudy , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 
    204203            ENDIF 
    205          END DO 
    206       ENDIF 
    207  
    208       !------------------------------------------- 
    209       ! Recover the properties from their contents 
    210       !------------------------------------------- 
    211       pato_i(:,:) = z0opw(:,:,1) * r1_e1e2t(:,:) * tmask(:,:,1) 
    212       DO jl = 1, jpl 
    213          pv_i (:,:,jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    214          pv_s (:,:,jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    215          psv_i(:,:,jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    216          poa_i(:,:,jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    217          pa_i (:,:,jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    218          DO jk = 1, nlay_s 
    219             pe_s(:,:,jk,jl) = z0es(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    220          END DO 
    221          DO jk = 1, nlay_i 
    222             pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    223          END DO 
    224          IF ( ln_pnd_H12 ) THEN 
    225             pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    226             pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     204            ! 
    227205         ENDIF 
     206 
     207         ! --- Recover the properties from their contents --- ! 
     208         DO jl = 1, jpl 
     209            pv_i (:,:,jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     210            pv_s (:,:,jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     211            psv_i(:,:,jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     212            poa_i(:,:,jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     213            pa_i (:,:,jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     214            DO jk = 1, nlay_s 
     215               pe_s(:,:,jk,jl) = z0es(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     216            END DO 
     217            DO jk = 1, nlay_i 
     218               pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     219            END DO 
     220            IF ( ln_pnd_H12 ) THEN 
     221               pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     222               pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     223            ENDIF 
     224         END DO 
     225         ! 
     226         ! derive open water from ice concentration 
     227         zati2(:,:) = SUM( pa_i(:,:,:), dim=3 ) 
     228         DO jj = 2, jpjm1 
     229            DO ji = fs_2, fs_jpim1 
     230               pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &                        !--- open water 
     231                  &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 
     232            END DO 
     233         END DO 
     234         CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T',  1. ) 
     235         ! 
     236         ! --- Ensure non-negative fields --- ! 
     237         !     Remove negative values (conservation is ensured) 
     238         !     (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 
     239         CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     240         ! 
     241         ! --- Ensure snow load is not too big --- ! 
     242         CALL Hsnow( zdt, pv_i, pv_s, pa_i, pa_ip, pe_s ) 
     243         ! 
    228244      END DO 
    229       ! 
    230       ! --- Ensure non-negative fields --- ! 
    231       !     Remove negative values (conservation is ensured) 
    232       !     (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 
    233       CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
    234       ! 
    235       ! --- Ensure snow load is not too big --- ! 
    236       CALL Hsnow( zdt, pv_i, pv_s, pa_i, pa_ip, pe_s ) 
    237245      ! 
    238246      IF( lrst_ice )   CALL adv_pra_rst( 'WRITE', kt )   !* write Prather fields in the restart file 
    … …  
    296304            DO ji = 1, jpi 
    297305               zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 
    298                zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt * e2u(ji,jj) / psm(ji,jj,jl) 
     306               zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl) 
    299307               zalfq        =  zalf * zalf 
    300308               zalf1        =  1.0 - zalf 
    … …  
    322330         DO jj = 2, jpjm1                      !  Flux from i+1 to i when u LT 0. 
    323331            DO ji = 1, fs_jpim1 
    324                zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt * e2u(ji,jj) / psm(ji+1,jj,jl)  
     332               zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl)  
    325333               zalg  (ji,jj) = zalf 
    326334               zalfq         = zalf * zalf 
    … …  
    465473            DO ji = fs_2, fs_jpim1 
    466474               zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) 
    467                zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt * e1v(ji,jj) / psm(ji,jj,jl) 
     475               zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl) 
    468476               zalfq        =  zalf * zalf 
    469477               zalf1        =  1.0 - zalf 
    … …  
    491499         DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0. 
    492500            DO ji = fs_2, fs_jpim1 
    493                zalf          = ( MAX(0._wp, -pvt(ji,jj) ) * pdt * e1v(ji,jj) ) / psm(ji,jj+1,jl)  
     501               zalf          = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl)  
    494502               zalg  (ji,jj) = zalf 
    495503               zalfq         = zalf * zalf 
    … …  
    645653      ! 
    646654      !                             !* allocate prather fields 
    647       ALLOCATE( sxopw(jpi,jpj,1)   , syopw(jpi,jpj,1)   , sxxopw(jpi,jpj,1)   , syyopw(jpi,jpj,1)   , sxyopw(jpi,jpj,1)   ,   & 
    648          &      sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) ,   & 
     655      ALLOCATE( sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) ,   & 
    649656         &      sxsn (jpi,jpj,jpl) , sysn (jpi,jpj,jpl) , sxxsn (jpi,jpj,jpl) , syysn (jpi,jpj,jpl) , sxysn (jpi,jpj,jpl) ,   & 
    650657         &      sxa  (jpi,jpj,jpl) , sya  (jpi,jpj,jpl) , sxxa  (jpi,jpj,jpl) , syya  (jpi,jpj,jpl) , sxya  (jpi,jpj,jpl) ,   & 
    … …  
    692699         !                                   !==========================! 
    693700         ! 
    694          IF( ln_rstart ) THEN   ;   id1 = iom_varid( numrir, 'sxopw' , ldstop = .FALSE. )    ! file exist: id1>0 
     701         IF( ln_rstart ) THEN   ;   id1 = iom_varid( numrir, 'sxice' , ldstop = .FALSE. )    ! file exist: id1>0 
    695702         ELSE                   ;   id1 = 0                                                  ! no restart: id1=0 
    696703         ENDIF 
    … …  
    728735            CALL iom_get( numrir, jpdom_autoglo, 'syyage', syyage ) 
    729736            CALL iom_get( numrir, jpdom_autoglo, 'sxyage', sxyage ) 
    730             !                                                        ! open water in sea ice 
    731             CALL iom_get( numrir, jpdom_autoglo, 'sxopw' , sxopw  ) 
    732             CALL iom_get( numrir, jpdom_autoglo, 'syopw' , syopw  ) 
    733             CALL iom_get( numrir, jpdom_autoglo, 'sxxopw', sxxopw ) 
    734             CALL iom_get( numrir, jpdom_autoglo, 'syyopw', syyopw ) 
    735             CALL iom_get( numrir, jpdom_autoglo, 'sxyopw', sxyopw ) 
    736737            !                                                        ! snow layers heat content 
    737738            DO jk = 1, nlay_s 
    … …  
    776777            sxsal = 0._wp   ;   sysal = 0._wp   ;   sxxsal = 0._wp   ;   syysal = 0._wp   ;   sxysal = 0._wp      ! ice salinity 
    777778            sxage = 0._wp   ;   syage = 0._wp   ;   sxxage = 0._wp   ;   syyage = 0._wp   ;   sxyage = 0._wp      ! ice age 
    778             sxopw = 0._wp   ;   syopw = 0._wp   ;   sxxopw = 0._wp   ;   syyopw = 0._wp   ;   sxyopw = 0._wp      ! open water in sea ice 
    779779            sxc0  = 0._wp   ;   syc0  = 0._wp   ;   sxxc0  = 0._wp   ;   syyc0  = 0._wp   ;   sxyc0  = 0._wp      ! snow layers heat content 
    780780            sxe   = 0._wp   ;   sye   = 0._wp   ;   sxxe   = 0._wp   ;   syye   = 0._wp   ;   sxye   = 0._wp      ! ice layers heat content 
    … …  
    825825         CALL iom_rstput( iter, nitrst, numriw, 'syyage', syyage ) 
    826826         CALL iom_rstput( iter, nitrst, numriw, 'sxyage', sxyage ) 
    827          !                                                           ! open water in sea ice 
    828          CALL iom_rstput( iter, nitrst, numriw, 'sxopw' , sxopw  ) 
    829          CALL iom_rstput( iter, nitrst, numriw, 'syopw' , syopw  ) 
    830          CALL iom_rstput( iter, nitrst, numriw, 'sxxopw', sxxopw ) 
    831          CALL iom_rstput( iter, nitrst, numriw, 'syyopw', syyopw ) 
    832          CALL iom_rstput( iter, nitrst, numriw, 'sxyopw', sxyopw ) 
    833827         !                                                           ! snow layers heat content 
    834828         DO jk = 1, nlay_s 
Note: See TracChangeset for help on using the changeset viewer.