Changeset 11732


Ignore:
Timestamp:
2019-10-18T18:10:30+02:00 (10 months ago)
Author:
clem
Message:

update trunk from #11730

Location:
NEMO/trunk/src/ICE
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/ICE/icedyn_adv_pra.F90

    r11632 r11732  
    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 
  • NEMO/trunk/src/ICE/icedyn_rdgrft.F90

    r11560 r11732  
    8686      !!                ***  ROUTINE ice_dyn_rdgrft_alloc *** 
    8787      !!------------------------------------------------------------------- 
    88       ALLOCATE( closing_net(jpij), opning(jpij)   , closing_gross(jpij),   & 
    89          &      apartf(jpij,0:jpl), hrmin(jpij,jpl), hraft(jpij,jpl)    , aridge(jpij,jpl),  & 
    90          &      hrmax(jpij,jpl), hi_hrdg(jpij,jpl)  , araft (jpij,jpl),  & 
     88      ALLOCATE( closing_net(jpij)  , opning(jpij)      , closing_gross(jpij) ,               & 
     89         &      apartf(jpij,0:jpl) , hrmin  (jpij,jpl) , hraft(jpij,jpl) , aridge(jpij,jpl), & 
     90         &      hrmax (jpij,jpl)   , hi_hrdg(jpij,jpl) , araft(jpij,jpl) ,                   & 
    9191         &      ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc ) 
    9292 
     
    137137      REAL(wp) ::   zfac                       ! local scalar 
    138138      INTEGER , DIMENSION(jpij) ::   iptidx        ! compute ridge/raft or not 
    139       REAL(wp), DIMENSION(jpij) ::   zdivu_adv     ! divu as implied by transport scheme  (1/s) 
    140139      REAL(wp), DIMENSION(jpij) ::   zdivu, zdelt  ! 1D divu_i & delta_i 
    141140      ! 
     
    175174         
    176175         ! just needed here 
    177          CALL tab_2d_1d( npti, nptidx(1:npti), zdivu   (1:npti)      , divu_i  ) 
    178176         CALL tab_2d_1d( npti, nptidx(1:npti), zdelt   (1:npti)      , delta_i ) 
    179177         ! needed here and in the iteration loop 
     178         CALL tab_2d_1d( npti, nptidx(1:npti), zdivu   (1:npti)      , divu_i) ! zdivu is used as a work array here (no change in divu_i) 
    180179         CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d  (1:npti,1:jpl), a_i   ) 
    181180         CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i   ) 
     
    187186            closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 
    188187            ! 
    189             ! divergence given by the advection scheme 
    190             !   (which may not be equal to divu as computed from the velocity field) 
    191             IF    ( ln_adv_Pra ) THEN 
    192                zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_rdtice 
    193             ELSEIF( ln_adv_UMx ) THEN 
    194                zdivu_adv(ji) = zdivu(ji) 
    195             ENDIF 
    196             ! 
    197             IF( zdivu_adv(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) )   ! make sure the closing rate is large enough 
    198             !                                                                                        ! to give asum = 1.0 after ridging 
     188            IF( zdivu(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu(ji) )   ! make sure the closing rate is large enough 
     189            !                                                                                ! to give asum = 1.0 after ridging 
    199190            ! Opening rate (non-negative) that will give asum = 1.0 after ridging. 
    200             opning(ji) = closing_net(ji) + zdivu_adv(ji) 
     191            opning(ji) = closing_net(ji) + zdivu(ji) 
    201192         END DO 
    202193         ! 
     
    215206               ato_i_1d   (ipti)   = ato_i_1d   (ji) 
    216207               closing_net(ipti)   = closing_net(ji) 
    217                zdivu_adv  (ipti)   = zdivu_adv  (ji) 
     208               zdivu      (ipti)   = zdivu      (ji) 
    218209               opning     (ipti)   = opning     (ji) 
    219210            ENDIF 
     
    259250               ELSE 
    260251                  iterate_ridging  = 1 
    261                   zdivu_adv  (ji) = zfac * r1_rdtice 
    262                   closing_net(ji) = MAX( 0._wp, -zdivu_adv(ji) ) 
    263                   opning     (ji) = MAX( 0._wp,  zdivu_adv(ji) ) 
     252                  zdivu      (ji) = zfac * r1_rdtice 
     253                  closing_net(ji) = MAX( 0._wp, -zdivu(ji) ) 
     254                  opning     (ji) = MAX( 0._wp,  zdivu(ji) ) 
    264255               ENDIF 
    265256            END DO 
     
    309300 
    310301      !                       ! Ice thickness needed for rafting 
    311       WHERE( pa_i(1:npti,:) > epsi20 )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:) 
     302      WHERE( pa_i(1:npti,:) > epsi10 )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:) 
    312303      ELSEWHERE                          ;   zhi(1:npti,:) = 0._wp 
    313304      END WHERE 
     
    328319      zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 ) 
    329320      ! 
    330       WHERE( zasum(1:npti) > epsi20 )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti) 
     321      WHERE( zasum(1:npti) > epsi10 )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti) 
    331322      ELSEWHERE                         ;   z1_asum(1:npti) = 0._wp 
    332323      END WHERE 
     
    454445      ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate.   
    455446      ! NOTE: 0 < aksum <= 1 
    456       WHERE( zaksum(1:npti) > epsi20 )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti) 
     447      WHERE( zaksum(1:npti) > epsi10 )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti) 
    457448      ELSEWHERE                          ;   closing_gross(1:npti) = 0._wp 
    458449      END WHERE 
     
    537528            IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN   ! only if ice is ridging 
    538529 
    539                IF( a_i_2d(ji,jl1) > epsi20 ) THEN   ;   z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) 
     530               IF( a_i_2d(ji,jl1) > epsi10 ) THEN   ;   z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) 
    540531               ELSE                                 ;   z1_ai(ji) = 0._wp 
    541532               ENDIF 
     
    595586               ! virtual salt flux to keep salinity constant 
    596587               IF( nn_icesal /= 2 )  THEN 
    597                   sirdg2(ji)     = sirdg2(ji)     - vsw * ( sss_1d(ji) - s_i_1d(ji) )        ! ridge salinity = s_i 
     588                  sirdg2(ji)     = sirdg2(ji)     - vsw * ( sss_1d(ji) - s_i_1d(ji) )       ! ridge salinity = s_i 
    598589                  sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoi * r1_rdtice  &  ! put back sss_m into the ocean 
    599590                     &                            - s_i_1d(ji) * vsw * rhoi * r1_rdtice     ! and get  s_i  from the ocean  
  • NEMO/trunk/src/ICE/iceitd.F90

    r11536 r11732  
    211211               CALL itd_glinear( zhb0(1:npti)  , zhb1(1:npti)  , h_ib_1d(1:npti)  , a_i_1d(1:npti)  ,  &   ! in 
    212212                  &              g0  (1:npti,1), g1  (1:npti,1), hL     (1:npti,1), hR    (1:npti,1)   )   ! out 
    213                   ! 
     213               ! 
    214214               ! Area lost due to melting of thin ice 
    215215               DO ji = 1, npti 
     
    218218                     ! 
    219219                     zdh0 =  h_i_1d(ji) - h_ib_1d(ji)                 
    220                      IF( zdh0 < 0.0 ) THEN      !remove area from category 1 
     220                     IF( zdh0 < 0.0 ) THEN      ! remove area from category 1 
    221221                        zdh0 = MIN( -zdh0, hi_max(1) ) 
    222222                        !Integrate g(1) from 0 to dh0 to estimate area melted 
     
    226226                           zx1    = zetamax 
    227227                           zx2    = 0.5 * zetamax * zetamax  
    228                            zda0   = g1(ji,1) * zx2 + g0(ji,1) * zx1                        ! ice area removed 
     228                           zda0   = g1(ji,1) * zx2 + g0(ji,1) * zx1                ! ice area removed 
    229229                           zdamax = a_i_1d(ji) * (1.0 - h_i_1d(ji) / h_ib_1d(ji) ) ! Constrain new thickness <= h_i                 
    230                            zda0   = MIN( zda0, zdamax )                                                  ! ice area lost due to melting  
    231                            !     of thin ice (zdamax > 0) 
     230                           zda0   = MIN( zda0, zdamax )                            ! ice area lost due to melting of thin ice (zdamax > 0) 
    232231                           ! Remove area, conserving volume 
    233232                           h_i_1d(ji) = h_i_1d(ji) * a_i_1d(ji) / ( a_i_1d(ji) - zda0 ) 
     
    349348      DO ji = 1, npti 
    350349         ! 
    351          IF( paice(ji) > epsi10  .AND. phice(ji) > 0._wp )  THEN 
     350         IF( paice(ji) > epsi10  .AND. phice(ji) > epsi10 )  THEN 
    352351            ! 
    353352            ! Initialize hL and hR 
  • NEMO/trunk/src/ICE/icevar.F90

    r11560 r11732  
    622622                  pv_s   (ji,jj,jl) = 0._wp 
    623623               ENDIF 
    624                IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 
     624               IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp .OR. pv_i(ji,jj,jl) <= 0._wp ) THEN 
    625625                  sfx_res(ji,jj)    = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt 
    626626                  psv_i  (ji,jj,jl) = 0._wp 
Note: See TracChangeset for help on using the changeset viewer.