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 13471 for NEMO/branches/2020/temporary_r4_trunk/src/ICE/iceistate.F90 – NEMO

Ignore:
Timestamp:
2020-09-16T11:40:06+02:00 (4 years ago)
Author:
smasson
Message:

r4_trunk: update iceistate.F90 for easy merge..., see #2523

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/temporary_r4_trunk/src/ICE/iceistate.F90

    r13470 r13471  
    170170      !------------------------------------------------------------------------ 
    171171      IF( ln_iceini ) THEN 
    172          !                             !---------------! 
    173          IF( nn_iceini_file == 1 )THEN ! Read a file   ! 
    174             !                          !---------------! 
    175             WHERE( ff_t(:,:) >= 0._wp )   ;   zswitch(:,:) = 1._wp 
    176             ELSEWHERE                     ;   zswitch(:,:) = 0._wp 
     172         ! 
     173         IF( Agrif_Root() ) THEN 
     174            !                             !---------------! 
     175            IF( nn_iceini_file == 1 )THEN ! Read a file   ! 
     176               !                          !---------------! 
     177               WHERE( ff_t(:,:) >= 0._wp )   ;   zswitch(:,:) = 1._wp 
     178               ELSEWHERE                     ;   zswitch(:,:) = 0._wp 
     179               END WHERE 
     180               ! 
     181               CALL fld_read( kt, 1, si ) ! input fields provided at the current time-step 
     182               ! 
     183               ! -- mandatory fields -- ! 
     184               zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1) * tmask(:,:,1) 
     185               zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1) * tmask(:,:,1) 
     186               zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1) * tmask(:,:,1) 
     187                
     188               ! -- optional fields -- ! 
     189               !    if fields do not exist then set them to the values present in the namelist (except for temperatures) 
     190               ! 
     191               ! ice salinity 
     192               IF( TRIM(si(jp_smi)%clrootname) == 'NOT USED' ) & 
     193                  &     si(jp_smi)%fnow(:,:,1) = ( rn_smi_ini_n * zswitch + rn_smi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 
     194               ! 
     195               ! temperatures 
     196               IF    ( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. & 
     197                  &    TRIM(si(jp_tms)%clrootname) == 'NOT USED' ) THEN 
     198                  si(jp_tmi)%fnow(:,:,1) = ( rn_tmi_ini_n * zswitch + rn_tmi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 
     199                  si(jp_tsu)%fnow(:,:,1) = ( rn_tsu_ini_n * zswitch + rn_tsu_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 
     200                  si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zswitch + rn_tms_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 
     201               ENDIF 
     202               IF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) & ! if T_s is read and not T_i, set T_i = (T_s + T_freeze)/2 
     203                  &     si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tms)%fnow(:,:,1) + 271.15 ) 
     204               IF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) & ! if T_su is read and not T_i, set T_i = (T_su + T_freeze)/2 
     205                  &     si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tsu)%fnow(:,:,1) + 271.15 ) 
     206               IF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) & ! if T_s is read and not T_su, set T_su = T_s 
     207                  &     si(jp_tsu)%fnow(:,:,1) = si(jp_tms)%fnow(:,:,1) 
     208               IF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) & ! if T_i is read and not T_su, set T_su = T_i 
     209                  &     si(jp_tsu)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1) 
     210               IF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) & ! if T_su is read and not T_s, set T_s = T_su 
     211                  &     si(jp_tms)%fnow(:,:,1) = si(jp_tsu)%fnow(:,:,1) 
     212               IF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) & ! if T_i is read and not T_s, set T_s = T_i 
     213                  &     si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1) 
     214               ! 
     215               ! pond concentration 
     216               IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) & 
     217                  &     si(jp_apd)%fnow(:,:,1) = ( rn_apd_ini_n * zswitch + rn_apd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) & ! rn_apd = pond fraction => rn_apnd * a_i = pond conc. 
     218                  &                              * si(jp_ati)%fnow(:,:,1)  
     219               ! 
     220               ! pond depth 
     221               IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) & 
     222                  &     si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 
     223               ! 
     224               ! pond lid depth 
     225               IF( TRIM(si(jp_hld)%clrootname) == 'NOT USED' ) & 
     226                  &     si(jp_hld)%fnow(:,:,1) = ( rn_hld_ini_n * zswitch + rn_hld_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 
     227               ! 
     228               zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1) * tmask(:,:,1) 
     229               ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1) * tmask(:,:,1) 
     230               zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1) * tmask(:,:,1) 
     231               ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1) * tmask(:,:,1) 
     232               zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1) * tmask(:,:,1) 
     233               zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1) * tmask(:,:,1) 
     234               zhlid_ini(:,:) = si(jp_hld)%fnow(:,:,1) * tmask(:,:,1) 
     235               ! 
     236               ! change the switch for the following 
     237               WHERE( zat_i_ini(:,:) > 0._wp )   ;   zswitch(:,:) = tmask(:,:,1)  
     238               ELSEWHERE                         ;   zswitch(:,:) = 0._wp 
     239               END WHERE 
     240               !                          !---------------! 
     241            ELSE                          ! Read namelist ! 
     242               !                          !---------------! 
     243               ! no ice if (sst - Tfreez) >= thresold 
     244               WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst )   ;   zswitch(:,:) = 0._wp  
     245               ELSEWHERE                                                                    ;   zswitch(:,:) = tmask(:,:,1) 
     246               END WHERE 
     247               ! 
     248               ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array 
     249               WHERE( ff_t(:,:) >= 0._wp ) 
     250                  zht_i_ini(:,:) = rn_hti_ini_n * zswitch(:,:) 
     251                  zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:) 
     252                  zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:) 
     253                  zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:) 
     254                  ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:) 
     255                  zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:) 
     256                  ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:) 
     257                  zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.  
     258                  zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:) 
     259                  zhlid_ini(:,:) = rn_hld_ini_n * zswitch(:,:) 
     260               ELSEWHERE 
     261                  zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:) 
     262                  zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:) 
     263                  zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:) 
     264                  zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:) 
     265                  ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:) 
     266                  zt_su_ini(:,:) = rn_tsu_ini_s * zswitch(:,:) 
     267                  ztm_s_ini(:,:) = rn_tms_ini_s * zswitch(:,:) 
     268                  zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 
     269                  zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:) 
     270                  zhlid_ini(:,:) = rn_hld_ini_s * zswitch(:,:) 
     271               END WHERE 
     272               ! 
     273            ENDIF 
     274             
     275            ! make sure ponds = 0 if no ponds scheme 
     276            IF ( .NOT.ln_pnd ) THEN 
     277               zapnd_ini(:,:) = 0._wp 
     278               zhpnd_ini(:,:) = 0._wp 
     279               zhlid_ini(:,:) = 0._wp 
     280            ENDIF 
     281             
     282            IF ( .NOT.ln_pnd_lids ) THEN 
     283               zhlid_ini(:,:) = 0._wp 
     284            ENDIF 
     285             
     286            !----------------! 
     287            ! 3) fill fields ! 
     288            !----------------! 
     289            ! select ice covered grid points 
     290            npti = 0 ; nptidx(:) = 0 
     291            DO_2D( 1, 1, 1, 1 ) 
     292               IF ( zht_i_ini(ji,jj) > 0._wp ) THEN 
     293                  npti         = npti  + 1 
     294                  nptidx(npti) = (jj - 1) * jpi + ji 
     295               ENDIF 
     296            END_2D 
     297             
     298            ! move to 1D arrays: (jpi,jpj) -> (jpi*jpj) 
     299            CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti)  , zht_i_ini ) 
     300            CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti)  , zht_s_ini ) 
     301            CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti)  , zat_i_ini ) 
     302            CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,1), ztm_i_ini ) 
     303            CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d (1:npti,1), ztm_s_ini ) 
     304            CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti)  , zt_su_ini ) 
     305            CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti)  , zsm_i_ini ) 
     306            CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti)  , zapnd_ini ) 
     307            CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d(1:npti)  , zhpnd_ini ) 
     308            CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d(1:npti)  , zhlid_ini ) 
     309             
     310            ! allocate temporary arrays 
     311            ALLOCATE( zhi_2d (npti,jpl), zhs_2d (npti,jpl), zai_2d (npti,jpl), & 
     312               &      zti_2d (npti,jpl), zts_2d (npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), & 
     313               &      zaip_2d(npti,jpl), zhip_2d(npti,jpl), zhil_2d(npti,jpl) ) 
     314 
     315            ! distribute 1-cat into jpl-cat: (jpi*jpj) -> (jpi*jpj,jpl) 
     316            CALL ice_var_itd( h_i_1d(1:npti)  , h_s_1d(1:npti)  , at_i_1d(1:npti),                  & 
     317               &              zhi_2d          , zhs_2d          , zai_2d         ,                  & 
     318               &              t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti),                  & 
     319               &              s_i_1d(1:npti)  , a_ip_1d(1:npti) , h_ip_1d(1:npti), h_il_1d(1:npti), & 
     320               &              zti_2d          , zts_2d          , ztsu_2d        ,                  & 
     321               &              zsi_2d          , zaip_2d         , zhip_2d        , zhil_2d ) 
     322 
     323            ! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl) 
     324            DO jl = 1, jpl 
     325               zti_3d(:,:,jl) = rt0 * tmask(:,:,1) 
     326               zts_3d(:,:,jl) = rt0 * tmask(:,:,1) 
     327            END DO 
     328            CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d   , h_i    ) 
     329            CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d   , h_s    ) 
     330            CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d   , a_i    ) 
     331            CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d   , zti_3d ) 
     332            CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d   , zts_3d ) 
     333            CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d  , t_su   ) 
     334            CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d   , s_i    ) 
     335            CALL tab_2d_3d( npti, nptidx(1:npti), zaip_2d  , a_ip   ) 
     336            CALL tab_2d_3d( npti, nptidx(1:npti), zhip_2d  , h_ip   ) 
     337            CALL tab_2d_3d( npti, nptidx(1:npti), zhil_2d  , h_il   ) 
     338 
     339            ! deallocate temporary arrays 
     340            DEALLOCATE( zhi_2d, zhs_2d, zai_2d , & 
     341               &        zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d, zhil_2d ) 
     342 
     343            ! calculate extensive and intensive variables 
     344            CALL ice_var_salprof ! for sz_i 
     345            DO jl = 1, jpl 
     346               DO_2D( 1, 1, 1, 1 ) 
     347                  v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl) 
     348                  v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl) 
     349                  sv_i(ji,jj,jl) = MIN( MAX( rn_simin , s_i(ji,jj,jl) ) , rn_simax ) * v_i(ji,jj,jl) 
     350               END_2D 
     351            END DO 
     352            ! 
     353            DO jl = 1, jpl 
     354               DO_3D( 1, 1, 1, 1, 1, nlay_s ) 
     355                  t_s(ji,jj,jk,jl) = zts_3d(ji,jj,jl) 
     356                  e_s(ji,jj,jk,jl) = zswitch(ji,jj) * v_s(ji,jj,jl) * r1_nlay_s * & 
     357                     &               rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus ) 
     358               END_3D 
     359            END DO 
     360            ! 
     361            DO jl = 1, jpl 
     362               DO_3D( 1, 1, 1, 1, 1, nlay_i ) 
     363                  t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl)  
     364                  ztmelts          = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K 
     365                  e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * & 
     366                     &               rhoi * (  rcpi  * ( ztmelts - t_i(ji,jj,jk,jl) ) + & 
     367                     &                         rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0), -epsi20 ) ) & 
     368                     &                       - rcp   * ( ztmelts - rt0 ) ) 
     369               END_3D 
     370            END DO 
     371             
     372            ! Melt ponds 
     373            WHERE( a_i > epsi10 )   ;   a_ip_eff(:,:,:) = a_ip(:,:,:) / a_i(:,:,:) 
     374            ELSEWHERE               ;   a_ip_eff(:,:,:) = 0._wp 
    177375            END WHERE 
     376            v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:) 
     377            v_il(:,:,:) = h_il(:,:,:) * a_ip(:,:,:) 
     378 
     379            ! specific temperatures for coupled runs 
     380            tn_ice(:,:,:) = t_su(:,:,:) 
     381            t1_ice(:,:,:) = t_i (:,:,1,:) 
    178382            ! 
    179             CALL fld_read( kt, 1, si ) ! input fields provided at the current time-step 
    180             ! 
    181             ! -- mandatory fields -- ! 
    182             zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1) * tmask(:,:,1) 
    183             zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1) * tmask(:,:,1) 
    184             zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1) * tmask(:,:,1) 
    185  
    186             ! -- optional fields -- ! 
    187             !    if fields do not exist then set them to the values present in the namelist (except for temperatures) 
    188             ! 
    189             ! ice salinity 
    190             IF( TRIM(si(jp_smi)%clrootname) == 'NOT USED' ) & 
    191                &     si(jp_smi)%fnow(:,:,1) = ( rn_smi_ini_n * zswitch + rn_smi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 
    192             ! 
    193             ! temperatures 
    194             IF    ( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. & 
    195                &    TRIM(si(jp_tms)%clrootname) == 'NOT USED' ) THEN 
    196                si(jp_tmi)%fnow(:,:,1) = ( rn_tmi_ini_n * zswitch + rn_tmi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 
    197                si(jp_tsu)%fnow(:,:,1) = ( rn_tsu_ini_n * zswitch + rn_tsu_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 
    198                si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zswitch + rn_tms_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 
    199             ENDIF 
    200             IF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) & ! if T_s is read and not T_i, set T_i = (T_s + T_freeze)/2 
    201                &     si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tms)%fnow(:,:,1) + 271.15 ) 
    202             IF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) & ! if T_su is read and not T_i, set T_i = (T_su + T_freeze)/2 
    203                &     si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tsu)%fnow(:,:,1) + 271.15 ) 
    204             IF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) & ! if T_s is read and not T_su, set T_su = T_s 
    205                &     si(jp_tsu)%fnow(:,:,1) = si(jp_tms)%fnow(:,:,1) 
    206             IF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) & ! if T_i is read and not T_su, set T_su = T_i 
    207                &     si(jp_tsu)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1) 
    208             IF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) & ! if T_su is read and not T_s, set T_s = T_su 
    209                &     si(jp_tms)%fnow(:,:,1) = si(jp_tsu)%fnow(:,:,1) 
    210             IF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) & ! if T_i is read and not T_s, set T_s = T_i 
    211                &     si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1) 
    212             ! 
    213             ! pond concentration 
    214             IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) & 
    215                &     si(jp_apd)%fnow(:,:,1) = ( rn_apd_ini_n * zswitch + rn_apd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) & ! rn_apd = pond fraction => rn_apnd * a_i = pond conc. 
    216                &                              * si(jp_ati)%fnow(:,:,1)  
    217             ! 
    218             ! pond depth 
    219             IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) & 
    220                &     si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 
    221             ! 
    222             ! pond lid depth 
    223             IF( TRIM(si(jp_hld)%clrootname) == 'NOT USED' ) & 
    224                &     si(jp_hld)%fnow(:,:,1) = ( rn_hld_ini_n * zswitch + rn_hld_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 
    225             ! 
    226             zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1) * tmask(:,:,1) 
    227             ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1) * tmask(:,:,1) 
    228             zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1) * tmask(:,:,1) 
    229             ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1) * tmask(:,:,1) 
    230             zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1) * tmask(:,:,1) 
    231             zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1) * tmask(:,:,1) 
    232             zhlid_ini(:,:) = si(jp_hld)%fnow(:,:,1) * tmask(:,:,1) 
    233             ! 
    234             ! change the switch for the following 
    235             WHERE( zat_i_ini(:,:) > 0._wp )   ;   zswitch(:,:) = tmask(:,:,1)  
    236             ELSEWHERE                         ;   zswitch(:,:) = 0._wp 
     383#if  defined key_agrif 
     384         ELSE 
     385 
     386            Agrif_SpecialValue    = -9999. 
     387            Agrif_UseSpecialValue = .TRUE. 
     388            CALL Agrif_init_variable(tra_iceini_id,procname=interp_tra_ice) 
     389            use_sign_north = .TRUE. 
     390            sign_north = -1. 
     391            CALL Agrif_init_variable(u_iceini_id  ,procname=interp_u_ice) 
     392            CALL Agrif_init_variable(v_iceini_id  ,procname=interp_v_ice) 
     393            Agrif_SpecialValue    = 0._wp 
     394            use_sign_north = .FALSE. 
     395            Agrif_UseSpecialValue = .FALSE. 
     396            ! lbc ????  
     397            ! Here we know : a_i, v_i, v_s, sv_i, oa_i, a_ip, v_ip, t_su, e_s, e_i 
     398            CALL ice_var_glo2eqv 
     399            CALL ice_var_zapsmall 
     400            CALL ice_var_agg(2) 
     401 
     402            ! Melt ponds 
     403            WHERE( a_i > epsi10 ) 
     404               a_ip_frac(:,:,:) = a_ip(:,:,:) / a_i(:,:,:) 
     405            ELSEWHERE 
     406               a_ip_frac(:,:,:) = 0._wp 
    237407            END WHERE 
    238             !                          !---------------! 
    239          ELSE                          ! Read namelist ! 
    240             !                          !---------------! 
    241             ! no ice if (sst - Tfreez) >= thresold 
    242             WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst )   ;   zswitch(:,:) = 0._wp  
    243             ELSEWHERE                                                                    ;   zswitch(:,:) = tmask(:,:,1) 
     408            WHERE( a_ip > 0._wp )       ! ???????     
     409               h_ip(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:) 
     410            ELSEWHERE 
     411               h_ip(:,:,:) = 0._wp 
    244412            END WHERE 
    245             ! 
    246             ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array 
    247             WHERE( ff_t(:,:) >= 0._wp ) 
    248                zht_i_ini(:,:) = rn_hti_ini_n * zswitch(:,:) 
    249                zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:) 
    250                zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:) 
    251                zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:) 
    252                ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:) 
    253                zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:) 
    254                ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:) 
    255                zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.  
    256                zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:) 
    257                zhlid_ini(:,:) = rn_hld_ini_n * zswitch(:,:) 
    258             ELSEWHERE 
    259                zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:) 
    260                zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:) 
    261                zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:) 
    262                zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:) 
    263                ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:) 
    264                zt_su_ini(:,:) = rn_tsu_ini_s * zswitch(:,:) 
    265                ztm_s_ini(:,:) = rn_tms_ini_s * zswitch(:,:) 
    266                zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 
    267                zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:) 
    268                zhlid_ini(:,:) = rn_hld_ini_s * zswitch(:,:) 
    269             END WHERE 
    270             ! 
    271          ENDIF 
    272  
    273          ! make sure ponds = 0 if no ponds scheme 
    274          IF ( .NOT.ln_pnd ) THEN 
    275             zapnd_ini(:,:) = 0._wp 
    276             zhpnd_ini(:,:) = 0._wp 
    277             zhlid_ini(:,:) = 0._wp 
    278          ENDIF 
    279  
    280          IF ( .NOT.ln_pnd_lids ) THEN 
    281             zhlid_ini(:,:) = 0._wp 
    282          ENDIF 
    283           
    284          !----------------! 
    285          ! 3) fill fields ! 
    286          !----------------! 
    287          ! select ice covered grid points 
    288          npti = 0 ; nptidx(:) = 0 
    289          DO_2D( 1, 1, 1, 1 ) 
    290             IF ( zht_i_ini(ji,jj) > 0._wp ) THEN 
    291                npti         = npti  + 1 
    292                nptidx(npti) = (jj - 1) * jpi + ji 
    293             ENDIF 
    294          END_2D 
    295  
    296          ! move to 1D arrays: (jpi,jpj) -> (jpi*jpj) 
    297          CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti)  , zht_i_ini ) 
    298          CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti)  , zht_s_ini ) 
    299          CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti)  , zat_i_ini ) 
    300          CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,1), ztm_i_ini ) 
    301          CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d (1:npti,1), ztm_s_ini ) 
    302          CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti)  , zt_su_ini ) 
    303          CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti)  , zsm_i_ini ) 
    304          CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti)  , zapnd_ini ) 
    305          CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d(1:npti)  , zhpnd_ini ) 
    306          CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d(1:npti)  , zhlid_ini ) 
    307  
    308          ! allocate temporary arrays 
    309          ALLOCATE( zhi_2d (npti,jpl), zhs_2d (npti,jpl), zai_2d (npti,jpl), & 
    310             &      zti_2d (npti,jpl), zts_2d (npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), & 
    311             &      zaip_2d(npti,jpl), zhip_2d(npti,jpl), zhil_2d(npti,jpl) ) 
    312           
    313          ! distribute 1-cat into jpl-cat: (jpi*jpj) -> (jpi*jpj,jpl) 
    314          CALL ice_var_itd( h_i_1d(1:npti)  , h_s_1d(1:npti)  , at_i_1d(1:npti),                  & 
    315             &              zhi_2d          , zhs_2d          , zai_2d         ,                  & 
    316             &              t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti),                  & 
    317             &              s_i_1d(1:npti)  , a_ip_1d(1:npti) , h_ip_1d(1:npti), h_il_1d(1:npti), & 
    318             &              zti_2d          , zts_2d          , ztsu_2d        ,                  & 
    319             &              zsi_2d          , zaip_2d         , zhip_2d        , zhil_2d ) 
    320  
    321          ! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl) 
    322          DO jl = 1, jpl 
    323             zti_3d(:,:,jl) = rt0 * tmask(:,:,1) 
    324             zts_3d(:,:,jl) = rt0 * tmask(:,:,1) 
    325          END DO 
    326          CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d   , h_i    ) 
    327          CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d   , h_s    ) 
    328          CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d   , a_i    ) 
    329          CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d   , zti_3d ) 
    330          CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d   , zts_3d ) 
    331          CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d  , t_su   ) 
    332          CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d   , s_i    ) 
    333          CALL tab_2d_3d( npti, nptidx(1:npti), zaip_2d  , a_ip   ) 
    334          CALL tab_2d_3d( npti, nptidx(1:npti), zhip_2d  , h_ip   ) 
    335          CALL tab_2d_3d( npti, nptidx(1:npti), zhil_2d  , h_il   ) 
    336  
    337          ! deallocate temporary arrays 
    338          DEALLOCATE( zhi_2d, zhs_2d, zai_2d , & 
    339             &        zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d, zhil_2d ) 
    340  
    341          ! calculate extensive and intensive variables 
    342          CALL ice_var_salprof ! for sz_i 
    343          DO jl = 1, jpl 
    344             DO_2D( 1, 1, 1, 1 ) 
    345                v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl) 
    346                v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl) 
    347                sv_i(ji,jj,jl) = MIN( MAX( rn_simin , s_i(ji,jj,jl) ) , rn_simax ) * v_i(ji,jj,jl) 
    348             END_2D 
    349          END DO 
    350          ! 
    351          DO jl = 1, jpl 
    352             DO_3D( 1, 1, 1, 1, 1, nlay_s ) 
    353                t_s(ji,jj,jk,jl) = zts_3d(ji,jj,jl) 
    354                e_s(ji,jj,jk,jl) = zswitch(ji,jj) * v_s(ji,jj,jl) * r1_nlay_s * & 
    355                   &               rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus ) 
    356             END_3D 
    357          END DO 
    358          ! 
    359          DO jl = 1, jpl 
    360             DO_3D( 1, 1, 1, 1, 1, nlay_i ) 
    361                t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl)  
    362                ztmelts          = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K 
    363                e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * & 
    364                   &               rhoi * (  rcpi  * ( ztmelts - t_i(ji,jj,jk,jl) ) + & 
    365                   &                         rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0), -epsi20 ) ) & 
    366                   &                       - rcp   * ( ztmelts - rt0 ) ) 
    367             END_3D 
    368          END DO 
    369  
    370          ! Melt ponds 
    371          WHERE( a_i > epsi10 )   ;   a_ip_eff(:,:,:) = a_ip(:,:,:) / a_i(:,:,:) 
    372          ELSEWHERE               ;   a_ip_eff(:,:,:) = 0._wp 
    373          END WHERE 
    374          v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:) 
    375          v_il(:,:,:) = h_il(:,:,:) * a_ip(:,:,:) 
    376            
    377          ! specific temperatures for coupled runs 
    378          tn_ice(:,:,:) = t_su(:,:,:) 
    379          t1_ice(:,:,:) = t_i (:,:,1,:) 
     413 
     414            tn_ice(:,:,:) = t_su(:,:,:) 
     415            t1_ice(:,:,:) = t_i (:,:,1,:) 
     416#endif 
     417         ENDIF ! Agrif_Root 
    380418         ! 
    381419         ! ice concentration should not exceed amax 
     
    440478      ENDIF 
    441479 
    442 !!clem: output of initial state should be written here but it is impossible because 
    443 !!      the ocean and ice are in the same file 
    444 !!      CALL dia_wri_state( 'output.init' ) 
     480      !!clem: output of initial state should be written here but it is impossible because 
     481      !!      the ocean and ice are in the same file 
     482      !!      CALL dia_wri_state( 'output.init' ) 
    445483      ! 
    446484   END SUBROUTINE ice_istate 
Note: See TracChangeset for help on using the changeset viewer.