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 12836 – NEMO

Changeset 12836


Ignore:
Timestamp:
2020-04-30T14:00:12+02:00 (4 years ago)
Author:
jcastill
Message:

Bug fix: when coupling with waves, the model crashed with an xios error saying that the mixed layer depth could not be written. This error is usually thrown when xios is trying to write NaNs?, but his was not the case. In this case the error looked to happen because xios did not finish writing a variable, but the variable was updated before xios finished. The do loop where this happened has been changed, and some cleaning of the branch done.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/UKMO/r12083_mix-lyr_diag/src/OCE/ZDF/zdfmxl.F90

    r12585 r12836  
    2828   PUBLIC   zdf_mxl       ! called by zdfphy.F90 
    2929 
    30    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld_tref  !: mixed layer depth at t-points - temperature criterion [m]  
    31    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   nmln       !: number of level in the mixed layer (used by TOP)  
    32    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld       !: mixing layer depth (turbocline)      [m]  
    33    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlp       !: mixed layer depth  (rho=rho0+zdcrit) [m]  
    34    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt      !: depth of the last T-point inside the mixed layer [m]  
    35    REAL(wp), PUBLIC, ALLOCATABLE,       DIMENSION(:,:) ::   hmld_zint  !: vertically-interpolated mixed layer depth   [m]    
    36    REAL(wp), PUBLIC, ALLOCATABLE,       DIMENSION(:,:) ::   htc_mld    ! Heat content of hmld_zint   
    37    LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)  ::   ll_found   ! Is T_b to be found by interpolation ?    
    38    LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)::   ll_belowml ! Flag points below mixed layer when ll_found=F 
     30   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hmld_tref  !: mixed layer depth at t-points - temperature criterion [m]  
     31   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   nmln       !: number of level in the mixed layer (used by TOP)  
     32   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hmld       !: mixing layer depth (turbocline)      [m]  
     33   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hmlp       !: mixed layer depth  (rho=rho0+zdcrit) [m]  
     34   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hmlpt      !: depth of the last T-point inside the mixed layer [m]  
     35   REAL(wp), PUBLIC, ALLOCATABLE,       DIMENSION(:,:,:) ::   hmld_zint  !: vertically-interpolated mixed layer depth   [m]    
     36   REAL(wp), PUBLIC, ALLOCATABLE,       DIMENSION(:,:,:) ::   htc_mld    ! Heat content of hmld_zint   
     37   LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    ::   ll_found   ! Is T_b to be found by interpolation ?    
     38   LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   ll_belowml ! Flag points below mixed layer when ll_found=F 
    3939 
    4040   REAL(wp), PUBLIC ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth 
    … …  
    5555   LOGICAL, PRIVATE :: mld_25h_write = .FALSE.  !Logical confirm 25h calculating/processing   
    5656   INTEGER, SAVE :: i_cnt_25h                   ! Counter for 25 hour means   
     57   INTEGER, PRIVATE :: nn_mld_diag = 0          ! number of diagnostics 
     58   INTEGER, PRIVATE, PARAMETER :: MAX_DIAG = 5  ! maximum number of diagnostics 
     59   LOGICAL, PRIVATE, DIMENSION(MAX_DIAG) :: cmld_zint, cmld_mld 
     60 
    5761   REAL(wp),SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: hmld_zint_25h 
    5862 
    … …  
    7074      zdf_mxl_alloc = 0      ! set to zero if no array to be allocated 
    7175      IF( .NOT. ALLOCATED( nmln ) ) THEN 
    72          ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), hmld_zint(jpi,jpj),       &  
    73                    htc_mld(jpi,jpj), &  
    74                    ll_found(jpi,jpj), ll_belowml(jpi,jpj,jpk), STAT= zdf_mxl_alloc ) 
     76         ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), hmld_zint(jpi,jpj, MAX_DIAG),       &  
     77                   htc_mld(jpi,jpj,MAX_DIAG), ll_found(jpi,jpj), ll_belowml(jpi,jpj,jpk), STAT= zdf_mxl_alloc ) 
    7578         ! 
    7679         ALLOCATE(hmld_tref(jpi,jpj)) 
    … …  
    190193      !!-----------------------------------------------------------------------------------  
    191194 
    192       TYPE(MXL_ZINT), INTENT(in)  :: sf 
     195      TYPE(MXL_ZINT), DIMENSION(MAX_DIAG), INTENT(in)  :: sf 
    193196 
    194197      ! Diagnostic criteria 
    … …  
    213216      REAL, DIMENSION(jpi,jpj)    :: zdelta_T      ! difference critereon  
    214217      REAL, DIMENSION(jpi,jpj)    :: zRHO1, zRHO2  ! Densities  
    215       INTEGER :: ji, jj, jk                             ! loop counter  
     218      INTEGER :: ji, jj, jk, jn                         ! loop counter  
    216219 
    217220      !!-------------------------------------------------------------------------------------  
    218221      !   
    219222      ! Unpack structure 
    220       nn_mld_type = sf%mld_type 
    221       rn_zref     = sf%zref 
    222       rn_dT_crit  = sf%dT_crit 
    223       rn_iso_frac = sf%iso_frac 
    224  
    225       ! Set the mixed layer depth criterion at each grid point  
    226       IF( nn_mld_type == 0 ) THEN 
    227          zdelta_T(:,:) = rn_dT_crit 
    228          zT(:,:,:) = rhop(:,:,:) 
    229       ELSE IF( nn_mld_type == 1 ) THEN 
    230          ppzdep(:,:)=0.0  
    231          call eos ( tsn(:,:,1,:), ppzdep(:,:), zRHO1(:,:) )  
     223      DO jn=1, nn_mld_diag 
     224         IF( cmld_zint(jn) .OR. cmld_mld(jn) ) THEN 
     225            nn_mld_type = sf(jn)%mld_type 
     226            rn_zref     = sf(jn)%zref 
     227            rn_dT_crit  = sf(jn)%dT_crit 
     228            rn_iso_frac = sf(jn)%iso_frac 
     229 
     230            ! Set the mixed layer depth criterion at each grid point  
     231            IF( nn_mld_type == 0 ) THEN 
     232               zdelta_T(:,:) = rn_dT_crit 
     233               zT(:,:,:) = rhop(:,:,:) 
     234            ELSE IF( nn_mld_type == 1 ) THEN 
     235               ppzdep(:,:)=0.0  
     236               call eos ( tsn(:,:,1,:), ppzdep(:,:), zRHO1(:,:) )  
    232237! Use zT temporarily as a copy of tsn with rn_dT_crit added to SST  
    233238! [assumes number of tracers less than number of vertical levels]  
    234          zT(:,:,1:jpts)=tsn(:,:,1,1:jpts)  
    235          zT(:,:,jp_tem)=zT(:,:,1)+rn_dT_crit  
    236          CALL eos( zT(:,:,1:jpts), ppzdep(:,:), zRHO2(:,:) )  
    237          zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0  
    238          ! RHO from eos (2d version) doesn't calculate north or east halo:  
    239          CALL lbc_lnk( 'zdfmxl', zdelta_T, 'T', 1. )  
    240          zT(:,:,:) = rhop(:,:,:)  
    241       ELSE  
    242          zdelta_T(:,:) = rn_dT_crit                       
    243          zT(:,:,:) = tsn(:,:,:,jp_tem)                            
    244       END IF  
    245  
    246       ! Calculate the gradient of zT and absolute difference for use later  
    247       DO jk = 1 ,jpk-2  
    248          zdTdz(:,:,jk)  =    ( zT(:,:,jk+1) - zT(:,:,jk) ) / e3w_n(:,:,jk+1)  
    249          zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) )  
    250       END DO  
    251  
    252       ! Find density/temperature at the reference level (Kara et al use 10m).           
    253       ! ik_ref is the index of the box centre immediately above or at the reference level  
    254       ! Find rn_zref in the array of model level depths and find the ref     
    255       ! density/temperature by linear interpolation.                                    
    256       DO jk = jpkm1, 2, -1  
    257          WHERE ( gdept_n(:,:,jk) > rn_zref )  
    258            ik_ref(:,:) = jk - 1  
    259            zT_ref(:,:) = zT(:,:,jk-1) + zdTdz(:,:,jk-1) * ( rn_zref - gdept_n(:,:,jk-1) )  
    260          END WHERE  
    261       END DO  
    262  
    263       ! If the first grid box centre is below the reference level then use the  
    264       ! top model level to get zT_ref  
    265       WHERE ( gdept_n(:,:,1) > rn_zref )   
    266          zT_ref = zT(:,:,1)  
    267          ik_ref = 1  
    268       END WHERE  
    269  
    270       ! The number of active tracer levels is 1 less than the number of active w levels  
    271       ikmt(:,:) = mbkt(:,:) - 1  
    272  
    273       ! Initialize / reset 
    274       ll_found(:,:) = .false. 
    275  
    276       IF ( rn_iso_frac - zepsilon > 0. ) THEN 
    277          ! Search for a uniform density/temperature region where adjacent levels           
    278          ! differ by less than rn_iso_frac * deltaT.                                       
    279          ! ik_iso is the index of the last level in the uniform layer   
    280          ! ll_found indicates whether the mixed layer depth can be found by interpolation  
    281          ik_iso(:,:)   = ik_ref(:,:)  
    282          DO jj = 1, nlcj  
    283             DO ji = 1, nlci  
     239               zT(:,:,1:jpts)=tsn(:,:,1,1:jpts)  
     240               zT(:,:,jp_tem)=zT(:,:,1)+rn_dT_crit  
     241               CALL eos( zT(:,:,1:jpts), ppzdep(:,:), zRHO2(:,:) )  
     242               zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0  
     243               ! RHO from eos (2d version) doesn't calculate north or east halo:  
     244               CALL lbc_lnk( 'zdfmxl', zdelta_T, 'T', 1. )  
     245               zT(:,:,:) = rhop(:,:,:)  
     246            ELSE  
     247               zdelta_T(:,:) = rn_dT_crit                       
     248               zT(:,:,:) = tsn(:,:,:,jp_tem)                            
     249            END IF  
     250 
     251            ! Calculate the gradient of zT and absolute difference for use later  
     252            DO jk = 1 ,jpk-2  
     253               zdTdz(:,:,jk)  =    ( zT(:,:,jk+1) - zT(:,:,jk) ) / e3w_n(:,:,jk+1)  
     254               zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) )  
     255            END DO  
     256 
     257            ! Find density/temperature at the reference level (Kara et al use 10m).           
     258            ! ik_ref is the index of the box centre immediately above or at the reference level  
     259            ! Find rn_zref in the array of model level depths and find the ref     
     260            ! density/temperature by linear interpolation.                                    
     261            DO jk = jpkm1, 2, -1  
     262               WHERE ( gdept_n(:,:,jk) > rn_zref )  
     263                 ik_ref(:,:) = jk - 1  
     264                 zT_ref(:,:) = zT(:,:,jk-1) + zdTdz(:,:,jk-1) * ( rn_zref - gdept_n(:,:,jk-1) )  
     265               END WHERE  
     266            END DO  
     267 
     268            ! If the first grid box centre is below the reference level then use the  
     269            ! top model level to get zT_ref  
     270            WHERE ( gdept_n(:,:,1) > rn_zref )   
     271               zT_ref = zT(:,:,1)  
     272               ik_ref = 1  
     273            END WHERE  
     274 
     275            ! The number of active tracer levels is 1 less than the number of active w levels  
     276            ikmt(:,:) = mbkt(:,:) - 1  
     277 
     278            ! Initialize / reset 
     279            ll_found(:,:) = .false. 
     280 
     281            IF ( rn_iso_frac - zepsilon > 0. ) THEN 
     282               ! Search for a uniform density/temperature region where adjacent levels           
     283               ! differ by less than rn_iso_frac * deltaT.                                       
     284               ! ik_iso is the index of the last level in the uniform layer   
     285               ! ll_found indicates whether the mixed layer depth can be found by interpolation  
     286               ik_iso(:,:)   = ik_ref(:,:)  
     287               DO jj = 1, nlcj  
     288                  DO ji = 1, nlci  
    284289!CDIR NOVECTOR  
    285                DO jk = ik_ref(ji,jj), ikmt(ji,jj)-1  
    286                   IF ( zmoddT(ji,jj,jk) > ( rn_iso_frac * zdelta_T(ji,jj) ) ) THEN  
    287                      ik_iso(ji,jj)   = jk  
    288                      ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) )  
    289                      EXIT  
    290                   END IF  
     290                     DO jk = ik_ref(ji,jj), ikmt(ji,jj)-1  
     291                        IF ( zmoddT(ji,jj,jk) > ( rn_iso_frac * zdelta_T(ji,jj) ) ) THEN  
     292                           ik_iso(ji,jj)   = jk  
     293                           ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) )  
     294                           EXIT  
     295                        END IF  
     296                     END DO  
     297                  END DO  
     298               END DO  
     299 
     300               ! Use linear interpolation to find depth of mixed layer base where possible  
     301               hmld_zint(:,:,jn) = rn_zref  
     302               DO jj = 1, jpj  
     303                  DO ji = 1, jpi  
     304                     IF (ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0) THEN  
     305                        zdz =  abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) )  
     306                        hmld_zint(ji,jj,jn) = gdept_n(ji,jj,ik_iso(ji,jj)) + zdz  
     307                     END IF  
     308                  END DO  
     309               END DO  
     310            END IF 
     311 
     312            ! If ll_found = .false. then calculate MLD using difference of zdelta_T     
     313            ! from the reference density/temperature  
     314 
     315! Prevent this section from working on land points  
     316            WHERE ( tmask(:,:,1) /= 1.0 )  
     317               ll_found = .true.  
     318            END WHERE  
     319 
     320            DO jk=1, jpk  
     321               ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= zdelta_T(:,:)   
     322            END DO  
     323 
     324! Set default value where interpolation cannot be used (ll_found=false)   
     325            DO jj = 1, jpj  
     326               DO ji = 1, jpi  
     327                  IF ( .not. ll_found(ji,jj) )  hmld_zint(ji,jj,jn) = gdept_n(ji,jj,ikmt(ji,jj))  
    291328               END DO  
    292329            END DO  
    293          END DO  
    294  
    295          ! Use linear interpolation to find depth of mixed layer base where possible  
    296          hmld_zint(:,:) = rn_zref  
    297          DO jj = 1, jpj  
    298             DO ji = 1, jpi  
    299                IF (ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0) THEN  
    300                   zdz =  abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) )  
    301                   hmld_zint(ji,jj) = gdept_n(ji,jj,ik_iso(ji,jj)) + zdz  
    302                END IF  
     330 
     331            DO jj = 1, jpj  
     332               DO ji = 1, jpi  
     333!CDIR NOVECTOR  
     334                  DO jk = ik_ref(ji,jj)+1, ikmt(ji,jj)  
     335                     IF ( ll_found(ji,jj) ) EXIT  
     336                     IF ( ll_belowml(ji,jj,jk) ) THEN                 
     337                        zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * SIGN(1.0, zdTdz(ji,jj,jk-1) )  
     338                        zdT  = zT_b - zT(ji,jj,jk-1)                                       
     339                        zdz  = zdT / zdTdz(ji,jj,jk-1)                                        
     340                        hmld_zint(ji,jj,jn) = gdept_n(ji,jj,jk-1) + zdz  
     341                        EXIT                                                    
     342                     END IF  
     343                  END DO  
     344               END DO  
    303345            END DO  
    304          END DO  
    305       END IF 
    306  
    307       ! If ll_found = .false. then calculate MLD using difference of zdelta_T     
    308       ! from the reference density/temperature  
    309   
    310 ! Prevent this section from working on land points  
    311       WHERE ( tmask(:,:,1) /= 1.0 )  
    312          ll_found = .true.  
    313       END WHERE  
    314   
    315       DO jk=1, jpk  
    316          ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= zdelta_T(:,:)   
    317       END DO  
    318   
    319 ! Set default value where interpolation cannot be used (ll_found=false)   
    320       DO jj = 1, jpj  
    321          DO ji = 1, jpi  
    322             IF ( .not. ll_found(ji,jj) )  hmld_zint(ji,jj) = gdept_n(ji,jj,ikmt(ji,jj))  
    323          END DO  
    324       END DO  
    325  
    326       DO jj = 1, jpj  
    327          DO ji = 1, jpi  
    328 !CDIR NOVECTOR  
    329             DO jk = ik_ref(ji,jj)+1, ikmt(ji,jj)  
    330                IF ( ll_found(ji,jj) ) EXIT  
    331                IF ( ll_belowml(ji,jj,jk) ) THEN                 
    332                   zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * SIGN(1.0, zdTdz(ji,jj,jk-1) )  
    333                   zdT  = zT_b - zT(ji,jj,jk-1)                                       
    334                   zdz  = zdT / zdTdz(ji,jj,jk-1)                                        
    335                   hmld_zint(ji,jj) = gdept_n(ji,jj,jk-1) + zdz  
    336                   EXIT                                                    
    337                END IF  
    338             END DO  
    339          END DO  
    340       END DO  
    341  
    342       hmld_zint(:,:) = hmld_zint(:,:)*tmask(:,:,1)  
     346 
     347            hmld_zint(:,:,jn) = hmld_zint(:,:,jn)*tmask(:,:,1)  
     348         END IF 
     349      END DO 
    343350      !   
    344351   END SUBROUTINE zdf_mxl_zint_mld 
    … …  
    355362      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    356363 
    357       INTEGER :: ji, jj, jk 
     364      INTEGER :: ji, jj, jk, jn 
    358365      INTEGER :: ikmax 
    359366      REAL(wp) :: zc, zcoef 
    … …  
    371378      ENDIF 
    372379 
    373       ! Find last whole model T level above the MLD 
    374       ilevel(:,:)   = 0 
    375       zthick_0(:,:) = 0._wp 
    376  
    377       DO jk = 1, jpkm1   
    378          DO jj = 1, jpj 
    379             DO ji = 1, jpi                     
    380                zthick_0(ji,jj) = zthick_0(ji,jj) + e3t_n(ji,jj,jk) 
    381                IF( zthick_0(ji,jj) < hmld_zint(ji,jj) )   ilevel(ji,jj) = jk 
     380      DO jn=1, nn_mld_diag 
     381         IF( cmld_mld(jn) ) THEN 
     382            ! Find last whole model T level above the MLD 
     383            ilevel(:,:)   = 0 
     384            zthick_0(:,:) = 0._wp 
     385 
     386            DO jk = 1, jpkm1   
     387               DO jj = 1, jpj 
     388                  DO ji = 1, jpi                     
     389                     zthick_0(ji,jj) = zthick_0(ji,jj) + e3t_n(ji,jj,jk) 
     390                     IF( zthick_0(ji,jj) < hmld_zint(ji,jj,jn) )   ilevel(ji,jj) = jk 
     391                  END DO 
     392               END DO 
     393               WRITE(numout,*) 'zthick_0(jk =',jk,') =',zthick_0(2,2) 
     394               WRITE(numout,*) 'gdepw_n(jk+1 =',jk+1,') =',gdepw_n(2,2,jk+1) 
    382395            END DO 
    383          END DO 
    384          WRITE(numout,*) 'zthick_0(jk =',jk,') =',zthick_0(2,2) 
    385          WRITE(numout,*) 'gdepw_n(jk+1 =',jk+1,') =',gdepw_n(2,2,jk+1) 
     396 
     397            ! Surface boundary condition 
     398            IF( ln_linssh ) THEN  ;   zthick(:,:) = sshn(:,:)   ;   htc_mld(:,:,jn) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1) 
     399            ELSE                  ;   zthick(:,:) = 0._wp       ;   htc_mld(:,:,jn) = 0._wp                                    
     400            ENDIF 
     401 
     402            ! Deepest whole T level above the MLD 
     403            ikmax = MIN( MAXVAL( ilevel(:,:) ), jpkm1 ) 
     404 
     405            ! Integration down to last whole model T level 
     406            DO jk = 1, ikmax 
     407               DO jj = 1, jpj 
     408                  DO ji = 1, jpi 
     409                     zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, ilevel(ji,jj) - jk + 1 ) , 1  )  )    ! 0 below ilevel 
     410                     zthick(ji,jj) = zthick(ji,jj) + zc 
     411                     htc_mld(ji,jj,jn) = htc_mld(ji,jj,jn) + zc * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk) 
     412                  END DO 
     413               END DO 
     414            END DO 
     415 
     416            ! Subsequent partial T level 
     417            zthick(:,:) = hmld_zint(:,:,jn) - zthick(:,:)   !   remaining thickness to reach MLD 
     418 
     419            DO jj = 1, jpj 
     420               DO ji = 1, jpi 
     421                  htc_mld(ji,jj,jn) = htc_mld(ji,jj,jn) + tsn(ji,jj,ilevel(ji,jj)+1,jp_tem)  &  
     422            &                      * MIN( e3t_n(ji,jj,ilevel(ji,jj)+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel(ji,jj)+1) 
     423               END DO 
     424            END DO 
     425 
     426            WRITE(numout,*) 'htc_mld(after) =',htc_mld(2,2,jn) 
     427 
     428            ! Convert to heat content 
     429            zcoef = rau0 * rcp 
     430            htc_mld(:,:,jn) = zcoef * htc_mld(:,:,jn) 
     431         END IF 
    386432      END DO 
    387  
    388       ! Surface boundary condition 
    389       IF( ln_linssh ) THEN  ;   zthick(:,:) = sshn(:,:)   ;   htc_mld(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)    
    390       ELSE                  ;   zthick(:,:) = 0._wp       ;   htc_mld(:,:) = 0._wp                                    
    391       ENDIF 
    392  
    393       ! Deepest whole T level above the MLD 
    394       ikmax = MIN( MAXVAL( ilevel(:,:) ), jpkm1 ) 
    395  
    396       ! Integration down to last whole model T level 
    397       DO jk = 1, ikmax 
    398          DO jj = 1, jpj 
    399             DO ji = 1, jpi 
    400                zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, ilevel(ji,jj) - jk + 1 ) , 1  )  )    ! 0 below ilevel 
    401                zthick(ji,jj) = zthick(ji,jj) + zc 
    402                htc_mld(ji,jj) = htc_mld(ji,jj) + zc * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk) 
    403             END DO 
    404          END DO 
    405       END DO 
    406  
    407       ! Subsequent partial T level 
    408       zthick(:,:) = hmld_zint(:,:) - zthick(:,:)   !   remaining thickness to reach MLD 
    409  
    410       DO jj = 1, jpj 
    411          DO ji = 1, jpi 
    412             htc_mld(ji,jj) = htc_mld(ji,jj) + tsn(ji,jj,ilevel(ji,jj)+1,jp_tem)  &  
    413       &                      * MIN( e3t_n(ji,jj,ilevel(ji,jj)+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel(ji,jj)+1) 
    414          END DO 
    415       END DO 
    416  
    417       WRITE(numout,*) 'htc_mld(after) =',htc_mld(2,2) 
    418  
    419       ! Convert to heat content 
    420       zcoef = rau0 * rcp 
    421       htc_mld(:,:) = zcoef * htc_mld(:,:) 
    422433 
    423434   END SUBROUTINE zdf_mxl_zint_htc 
    … …  
    437448      INTEGER :: jn 
    438449 
    439       INTEGER :: nn_mld_diag = 0    ! number of diagnostics 
    440  
    441450      CHARACTER(len=1) :: cmld 
    442  
    443       TYPE(MXL_ZINT) :: sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5 
    444       TYPE(MXL_ZINT), SAVE, DIMENSION(5) ::   mld_diags 
    445  
    446       NAMELIST/namzdf_mldzint/ nn_mld_diag, sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5 
     451      TYPE(MXL_ZINT), SAVE, DIMENSION(MAX_DIAG) ::   mld_diags 
     452 
     453      NAMELIST/namzdf_mldzint/ nn_mld_diag, mld_diags 
    447454 
    448455      !!---------------------------------------------------------------------- 
    … …  
    458465         IF(lwm) WRITE ( numond, namzdf_mldzint ) 
    459466 
    460          IF( nn_mld_diag > 5 )   CALL ctl_stop( 'STOP', 'zdf_mxl_ini: Specify no more than 5 MLD definitions' ) 
    461  
    462          mld_diags(1) = sn_mld1 
    463          mld_diags(2) = sn_mld2 
    464          mld_diags(3) = sn_mld3 
    465          mld_diags(4) = sn_mld4 
    466          mld_diags(5) = sn_mld5 
    467  
    468          IF( lwp .AND. (nn_mld_diag > 0) ) THEN 
    469             WRITE(numout,*) '=============== Vertically-interpolated mixed layer ================' 
    470             WRITE(numout,*) '(Diagnostic number, nn_mld_type, rn_zref, rn_dT_crit, rn_iso_frac)' 
     467         WRITE(cmld,'(I1)') MAX_DIAG  
     468         IF( nn_mld_diag > MAX_DIAG )   CALL ctl_stop( 'STOP', 'zdf_mxl_ini: Specify no more than ', 'cmld', ' MLD definitions' ) 
     469 
     470         cmld_zint=.false. 
     471         cmld_mld=.false. 
     472         IF( nn_mld_diag > 0 ) THEN 
     473            IF( lwp ) THEN 
     474               WRITE(numout,*) '=============== Vertically-interpolated mixed layer ================' 
     475               WRITE(numout,*) '(Diagnostic number, nn_mld_type, rn_zref, rn_dT_crit, rn_iso_frac)' 
     476            END IF 
     477 
    471478            DO jn = 1, nn_mld_diag 
    472                WRITE(numout,*) 'MLD criterion',jn,':' 
    473                WRITE(numout,*) '    nn_mld_type =', mld_diags(jn)%mld_type 
    474                WRITE(numout,*) '    rn_zref ='    , mld_diags(jn)%zref 
    475                WRITE(numout,*) '    rn_dT_crit =' , mld_diags(jn)%dT_crit 
    476                WRITE(numout,*) '    rn_iso_frac =', mld_diags(jn)%iso_frac 
     479               ! Check if the diagnostics is being written to the output 
     480               WRITE(cmld,'(I1)') jn 
     481               IF( iom_use( "mldzint_"//cmld ) ) cmld_zint(jn)=.true.  
     482               IF( iom_use( "mldhtc_"//cmld ) )  cmld_mld(jn) =.true.  
     483 
     484               IF( lwp ) THEN 
     485                  WRITE(numout,*) 'MLD criterion',jn,':' 
     486                  WRITE(numout,*) '    nn_mld_type =', mld_diags(jn)%mld_type 
     487                  WRITE(numout,*) '    rn_zref ='    , mld_diags(jn)%zref 
     488                  WRITE(numout,*) '    rn_dT_crit =' , mld_diags(jn)%dT_crit 
     489                  WRITE(numout,*) '    rn_iso_frac =', mld_diags(jn)%iso_frac 
     490               END IF 
    477491            END DO 
    478492            WRITE(numout,*) '====================================================================' 
    … …  
    481495 
    482496      IF( nn_mld_diag > 0 ) THEN 
     497         CALL zdf_mxl_zint_mld( mld_diags ) 
     498         CALL zdf_mxl_zint_htc( kt ) 
     499 
    483500         DO jn = 1, nn_mld_diag 
    484501            WRITE(cmld,'(I1)') jn 
    485             IF( iom_use( "mldzint_"//cmld ) .OR. iom_use( "mldhtc_"//cmld ) ) THEN 
    486                CALL zdf_mxl_zint_mld( mld_diags(jn) ) 
    487  
    488                IF( iom_use( "mldzint_"//cmld ) ) THEN 
    489                   CALL iom_put( "mldzint_"//cmld, hmld_zint(:,:) ) 
    490                ENDIF 
    491  
    492                IF( iom_use( "mldhtc_"//cmld ) )  THEN 
    493                   CALL zdf_mxl_zint_htc( kt ) 
    494                   CALL iom_put( "mldhtc_"//cmld , htc_mld(:,:)   ) 
    495                ENDIF 
     502            IF( cmld_zint(jn) ) THEN 
     503               CALL iom_put( "mldzint_"//cmld, hmld_zint(:,:,jn) ) 
     504            ENDIF 
     505 
     506            IF( cmld_mld(jn) )  THEN 
     507               CALL iom_put( "mldhtc_"//cmld , htc_mld(:,:,jn)   ) 
    496508            ENDIF 
    497509         END DO 
Note: See TracChangeset for help on using the changeset viewer.